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1. INTRODUCTION 

The purpose of this report is to discover and quantify the systematic 
errors in the algorithms employed by NUWES in their short base line 
underwater position location system. Systematic error can have a number of 
sources and previous works [6,7] have treated other issues. The present work 
deals with the algorithms used in sound ray tracing. There are a number of 
aspects to be treated, and some background is necessary in order to explain 
them. 

Figure 1 contains a schematic diagram of a short base line hydrophonic 
array and of the signals it may receive. Sharply pulsed signals, or pings, are 
sent by the sound source vehicles (surface craft, submarine, torpedo) and they 
are received by the four transducers (called the X, Y, Z and C-phones) of the 
array. These four hydrophones form a right angled coordinate system with 
origins at the C-phone and arm lengths D (=30 feet) to each of the other three 
phones. 

The ray paths from a source to the four receivers are synchronously timed 

_7 

with great precision (10 secs) and the differentials of arrival times are used 
to construct the direction of the source. But due to variability of the speed of 
sound at various water depths, the ray paths themselves are not straight lines. 
Also the paths may change from day to day as the water depth-velocity profile 
changes. Knowledge of the current profile allows one to perform a ray tracing 
computation. Recovery of the three-dimensional position of the sound 
source is accomplished by reconstructing the ray path and following it for the 
given amount of time. Each source vehicle has a phase coded "ping" so that 
its signals can be discriminated from those of other vehicles. 
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Each hydrophonic array is placed on the sea floor by lowering it over the 
side of a utility ship. Each has a self-leveling capability. When finally resting 
on the bottom, they are not perfectly level, but the X and Y arms have tilt 
meters that measure the angles that are made with the horizontal. An array 
surveying activity checks these angles and measures the rotation of the 
vertical. Thus the local coordinate system can be reconciled with the master 
range coordinate system. 

The currently employed ray tracing methodology partitions the water 
column into a number of equal thickness layers and treats the speed of sound 
as constant in each layer. This leads to the use of isospeed ray tracing [1, 5]. 
Presently the layer thicknesses at Nanoose are 25 feet. In those instances for 
which the measured depth-velocity profile does not go as deep as the array, a 
constant speed extrapolation is used. In effect the thickness of the deepest 
layer is larger, perhaps 50 or 75 feet. 

Now the issues can be detailed: 

i) How accurate is iso-speed raytracing using a 25 foot water layer 
increment? 

ii) What is the effect when it is necessary to use a thicker deepest 
layer? 

iii) How well are the ray tracing directions and transit times 
determined? 

iv) What effect does the various depth velocity profiles have upon 
the answers to i), ii) and iii)? 

The treatment of these questions requires a valid sound ray construction 
methodology and a representative set of depth-velocity profiles. To satisfy the 
latter requirement we have selected twelve experimental days at the Nanoose 
range spanning the period May 1988 to June 1989. They are presented (in 
graphical form) in Appendix A. We note that the more interesting ones seem 
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to appear in the Spring. Two of the experimental days, 23-24 April, 1989 are 
consecutive. This allows us a glimpse into the question of day to day 
variability. 

To treat the former requirement, we have developed an isogradient ray 
fitting algorithm. It fits a sound ray connecting two given points (in the 
horizontal-vertical plane) assuming direct path propagation. The depth 
velocity profile is partitioned into five foot equi thickness layers and within 
each layer the slope of sound speed vs. depth is constant. Thus the DV profile 
is represented as a continuous function consisting of a sequence of straight 
line segments. Within each layer the ray path is a circle arc because of Snell's 
law, [2,3,8]. The outputs of this algorithm are the angles that the ray makes 
with the horizontal at each of the endpoints, and the transit time of the ray 
from the initial point to the final one. 

Section 2 of the report presents some theoretical material and formulas. 
The distinctions between isospeed and isogradient ray tracing are explained. 
Snell's law is introduced and supported. 

Section 3 of the report deals with the accuracy of pure ray tracing in two 
dimensions, the horizontal-vertical plane. Computations are made for a 
number of initial angle-transit time pairs and for all twelve depth velocity 
profiles. Errors from this source are generally small but can be as large as a 
foot or more. The situation is more difficult when extrapolation of the water 
column is necessary. This occurs when the sounding does not extend as deep 
as the receiver. In such cases we cannot be definite about the nature of the 
errors, but several equally defensible methods lead to results that disagree by 
five or ten feet and even more. We state that there is a problem with 
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extrapolation. The soundings should be made at the deepest part of the range 
and to full depth. Failing that, a careful development of an extrapolation 
policy should be made. 

Section 4 deals with the three dimensional problem of locating the 
position of the source based upon the transit times to each of the four 
hydrophones. Hence the question is one of finding the azimuth and initial 
elevation angles of the ray, and a matching transit time to stop the ray tracing 
algorithm. Also there are confounding sources of variability. The depth- 
velocity profile plays a role, as mentioned before. But the arrays themselves 
have directional properties that would interact with the algorithm even if 
they were fully level and aligned with the range. The fact that the arrays are 
tilted and rotated in a variety of ways has contributed to the puzzle of 
interpreting the mismatches in the array overlap areas. Taken altogether it is 
shown that systematic error from these sources can be as large as ten or 
twelve feet. Moreover errors of this magnitude are unnecessary. An 
alternative method is proposed which can reduce them by at least an order of 
magnitude. 

The conclusions are summarized in Section 5. Section 6, an addendum, 
addresses an issue raised in the review process. Also a number of appendices 
are included. They hold the depth-velocity profiles, supporting mathematical 
details, details of the algorithms, and the source code for the FORTRAN 
programs. 

A brief general statement of conclusions is as follows: 

i) The error in iso-speed raytracing is an increasing function of 
horizontal range, but is seldom more than one foot. 
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ii) The error due to constant speed extrapolation in the deepest 
layer can range up to 10 or more feet. 

iii) The error due to initializing the ray tracing is a periodic 
function of the azimuth direction and can be substantial for 
tilted arrays and at the greater horizontal ranges. The effect of 
the determination of azimuth is especially noticeable. 

Greater details are presented as the various issues are developed. 

2. PERTINENT ITEMS FROM RAY TRACING 

The sound source pings and sends out an isotropic wave front, which is a 
fixed phase point on the pressure cycle. A receiving transducer measures the 
time of arrival of the wave front. A ray is a path normal to the wave front 
and extending in space back to the source. Our goal is to trace a ray from the 
receiver for a fixed amount of time and thereby locate the source. To do this 
we must construct the azimuth and elevation angles of the ray at the receiver 
and then bend the ray back through the various speed layers until the 
measured transit time is consumed. 

If the speed of sound in water were constant then the ray path would be a 
straight line. But it is not. Speed is a function of temperature, pressure and 
salinity. These variables interact in interesting ways for the water layer that 
affects our problem. Speed is not necessarily a monotone function of depth. 
Conditions change with time, and water sounding drops are made daily. 
They provide a depth- velocity profile which is assumed to be fixed for the 
entire day's exercises. Further, these values are assumed constant throughout 
each horizontal plane; i.e., the field is homogeneous. 

Our immediate goal is to justify the use of a ray invariant in a 
horizontal-vertical plane and to establish the circular arc nature of ray paths 
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in water layers for which the speed of sound is a straight line function of 
depth. 




s 2 

si 



Let us begin with 
Snell's law. Consider two 
adjacent layers with speed 
Si in the lower layer and 
speed S2 in the upper. The 

ray enters the lower layer at elevation angle 0i and the upper layer with angle 
02- Given si and S2 let us find the relationship between 0] and 02 that will 
minimize the transit time from (0,0) to (c,a+b). 

Proposition. For a ray to traverse from (0,0) to (c,a+b) in minimum time, 
we must have the relationship 

cos(0i) cos(02) 



S2 



( 2 . 1 ) 



Proof. The transit time of a path from (0,0) to (c,a+b) that goes through 
(x,a) is given by 
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T(x) = — Va 2 + x 2 + — Vb 2 + (c-x) 2 . 

Si S2 



Further, it has derivative 

lx 1 c-x 

T (x) = _ . 

s i Va 2 + x 2 S2 Vb 2 + (c-x) 2 

The relationship (2.1) is a consequent of setting T'(x) = 0. 

Now suppose the point c is not fixed but variable. It follows that a ray 
entering the lower layer at an angle 0i (with the horizontal) will seek the path 
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of minimum transit time and exit the upper layer at an angle of 62 ; the two 
angles are related by ( 2 . 1 ). 

Next suppose that a number of layers are stacked vertically; within each 
layer the speed of sound is constant. The relationship ( 2 . 1 ) must hold for 
every successive pair of layers and hence the ratio 

cos( 0 ) 



rv = 



( 2 . 2 ) 



must be constant for the entire ray path. This value characterizes the ray path 
and is called the ray invariant. This is Snell's law. 

Consider the vertical plane containing the source and the receiver. Call 
this is (h,z) plane with the depth z taken as positive downward. Our ray 
tracing problem is two dimensional in this plane. Given the depth velocity 
profile we need only the elevation angle at the receiver and the transit time 
to locate the source. 

The depth- velocity profile can be approximated with a series of thin layers 
each having constant (internal) sound speed. Thus ray tracing can be enacted 
using such an approximation. This approach is called isospeed ray tracing [ 1 ], 

A sharper approximation is available through the use of water layers 
whose sound speed structure can be represented with a linear function of 
depth 

u(z) = oo + t)]Z (2.3) 
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Proposition. If (2.3) holds then 
the ray path is a circle arc of 
radius z m + \)o/"Ui. If Uj > 0 then 
the arc is a distance Oo/Oi above 
z = 0. 

Proof. Referring to the diagram, 
let (h,z) be an arbitrary point on 
the ray path and 0 is the angle of the tangent line. 

Because of Snell's law we must have 

cos(9) 1 

UO+UiZ = -O0+ -OlZmax = ™ 

and the ray path can be described parametrically in terms of 

{h(0), z(0)} for 0 < 0 

The radius of curvation can be found from the general formula 

3/2 



r = ■ 



([h'(0)P + [z'(0)]2) 



dh 



h'(0)z"(0) - h"(0)z'(0) | ‘ 

Using (2.3), = cot(0), and implicit differentiation we find 

z'(0) = -sin(0)/i)irv h'(0) = -cos(0)/\)irv 

z"(0) = -cos(0)/uirv h"(0) = sin(0)/oirv 

and hence 

r = 1/uirv =— + z max 



(2.4) 
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which is a constant, independent of 0. Hence the path is a circle arc; the 
radius is as specified; and the center of the circle is on a line above zero at a 
distance uo/ui- This proof follows those found in [4,8]. 

If \)i < 0, then the [h(0), z(0)] curve is concave, still a circle arc, and the 
circle's center is still on the line z = -\)o/t)i. But now this line is below z=0. If 
\)i = 0, the circle radius is infinite, the sound speed is constant, and the ray 
path is a straight line. 

Now, the numerical construction of the ray path can also be accomplished 
by representing each layer's sound speed as a straight line segment (function 
of depth) and piecing together the consequent circle arcs. Since each layer has 
a constant speed gradient with depth, this is called isogradient ray tracing, [1]. 

The question of efficiency of the two methods, isospeed and isogradient 
ray tracing, is really a question of how well the depth-speed function in a layer 
can be represented by a constant on the one hand, or by a linear function on 
the other. For thick layers there may be oscillations that make the choice 
difficult. For thin layers it seems that the straight line fit should perform 
better. 

The algorithm for isogradient ray tracing is presented in Appendix C. A 
corresponding algorithm for isospeed ray tracing can be extracted from it by 
making a number of deletions. Fortran source codes for each are shown in 
Appendix G. Inputs for these algorithms include the depth-velocity table; the 
layer depths; the depth of the receiver; the elevation (layer entrance) angle at 
the receiver; the ray transit time. The outputs are the horizontal and vertical 
end points of the ray, and the final elevation (exit layer) angle of the ray. This 
latter quantity is needed in the timing synchronization model, [7], 
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For isospeed ray tracing the pertinent formulae are 



Ah = Az cot(0) 

At = Az/ X) sin(0) 

where 0 is the angle that the ray enters the layer; Az is the layer thickness; Ah 
is the horizontal distance transversed in the layer; At is the transit time 
through the layer. The next layer entrance angle is computed from the ray 
invariant equation (2.1). 

For isogradient tracing the pertinent formulas are more complicated. For 
a ray that enters a layer at depth zo; horizontal displacement ho; and angle 0o 
we must first compute the coordinate (qi, qo) of the center of the circle arc 
Coi*0) 

q 2 = -uo/vi 

(2.5) 

qi = h 0 + (q 2 -zo) sin(0o)/cos(0o) 
and the radius of the arc 

r = signum (q 2 ) (q 2 - zo)/cos(0o). (2.6) 

The new horizontal displacement is 

hi = qi - signum(q 2 ) r sin(0 o ) (2.7) 



r as 

and the increase in transit time is the line integral At = — along the 

J uq + Ms) 

circular path. This integral is most easily managed using ds = -\J(dz)~ +(dh)^ = 
dd / i )\rv and Uq + t>iz = cos(0)/rv so that 



At 



”1 



t d6 - 1 {in 


l + sin(0]) 


-In 


l + sin(6 0 )l] 


J o cos(0) ■ V, | 


cos(6i) 


c °s(6 0 ) Jj 
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and the layer exit angle is computed from the ray invariant equation .2.1), 

cos(0i) = rvu(zi) 



The layer exit angle is the entrance angle for the next layer. 

The above equations form the heart of direct path ray tracing. The 
organizational questions that arise when developing a ray tracing algorithm 
are treated in Appendix C. Basically one proceeds upwards through the layers 
until the specified total transit time is consumed. An end correction is 
normally necessary because of the requirement to stop part way through a 
layer. 

3. QUALITY OF ISOSPEED RAY TRACING 

We are concerned with the quality of the currently employed isospeed ray 
tracing algorithm, which uses a uniform water layer thickness of 25 feet. The 
receiver depths range from about 1100 to 1350 feet. The elevation angle can 
range from 90° (directly overhead) to some rather small but positive values. 
Of course a variety of water columns (depth-velocity profiles) can be 
encountered. We have selected twelve (see Appendix A) spanning the period 
May, 1988 to July 1989. 

Our first problem is to establish a standard ray to serve as a basis for 
comparison. To this end we are limited by the resolution of the water 
column data available. The information depicted graphically in Appendix A 
provides sound speed averages for every five feet. That is, at level / the 
corresponding velocity value represents 






^+5 




(3.1) 
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and we have no information concerning the amount of variability that may 
exist within the layer. It is presumed small and is neglected. (A model for 
assessing such variability is presented in Appendix D, and this author is 
concerned about the issue for small entrance angles.) 

The most expedient standard available is to employ the ray established by 
isogradient ray tracing utilizing the five foot layer thicknesses with the 
straight line segments as depicted (and exaggerated) in Figure 2. These rays 
are used to judge the rays formed by the isospeed method with 25 foot layer 
thickness. The corresponding depth-velocity table is formed by partitioning 
the {u/} into consecutive sets of size 5 and, within each set, average the five 
values. 

With the above as background, the remainder of this section deals with 
numerical comparisons treating three issues: the computational noise 

generated by the processing of a large number of layers on the computer; the 
basic precision of the current isospeed ray tracing; the effects of extrapolation 
policies when the measured water column does not go sufficiently deep. 

i. The adopted standard generally processes over 200 water layers (20 for 
each 100 feet separating source and receiver). The buildup of computational 
noise can be checked by using an artificial but linear depth velocity profile: 
An exact ray can be developed from the theory presented in the previous 
section and applied to a single layer of great thickness. Then the isogradient 
programs can attempt to match this ray by tracing through the usual five foot 
layers. This was done for the depth layer 150, 1200 ft. and intercepts 4250, 4800, 
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Depth 




Figure 2. Schematic Diagram Comparing Isogradient and Isospeed 
Representations of Depth-Velocity Information 
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4850 ft/sec matched with slopes .1, .05, .01 ft/secft respectively. The two 
methods produce virtually identical results. All computations are in double 
precision arithmetic. 

ii. The ray fitting technique produces transit time, and entrance and exit 
angles of the ray connecting two points (source and receiver coordinates) in 
the horizontal-depth plane. When isogradient ray tracing is applied starting 
at the receiver, using the entrance angle and stopping when the transit term 
is consumed, then the coordinates of the source are reproduced to within 
some small preassigned value e (we used e = 10 -6 ). 

The accuracy of isospeed raytracing using 25 foot layer increments is 
compared against isogradient ray tracing using five foot-layer increments. A 
fixed set of entrance angle (9 0 in radians) and transit time (to in seconds) pairs 
have been selected. Generally they produce horizontal distances of 1000 to 
6000 feet and depths of less than 50 to over 500 feet. The distance, d, between 
the two versions of sound source is compared for each input pair (Go, to) and 
each of the twelve DV profiles. The values are recorded in Table 1 along with 
the source coordinates (h c , z c ) for currently utilized isospeed methodology, 
and (h s , z s ) for the standard (isogradient). An additional computation (h g/ z g ) 
was performed using isogradient ray tracing with 25 ft layer thicknesses. (The 
purpose was to provide an indication of the relative importance of layer 
thickness and the two types of ray tracing.) In all cases the receiver depth is 
1300 feet. 

The errors are computed using 




(3.2) 
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with the subscript r replaced by c and g, respectively, in order to identify the 
current and 25 ft isogradient computations. 

Examination of Table 1 shows that the errors, d, grow with range and the 
majority of the total error is in the vertical component. The effect of varying 
DV profiles is quite noticeable. The larger errors can exceed one foot. 
Generally isogradient ray tracing with 25-foot layer thicknesses is noticeably 
better than the current methods. 

iii. Some measure of the effect of extrapolation techniques is given in 
Table 2. The same cases are developed as found in Table 1, but the sensor 
depth has increased to 1350 feet. The DV profiles seldom goes below 1300 feet 
so extrapolation of speed information is necessary. There are quite a few 
arrays deeper than 1300 feet and several (2, 3, 7, 13, 14 see Table B-l) 
considerably so. Moreover, the C-phones are even deeper, see Table B-2. 
(Fourteen of the C-phones are deeper than 1325 ft.) The methods of 
extrapolation are explained in Appendix A along with the visual effect of 
their use. In some instances the visual effect is appealing and in others it is 
not. See the insets in Appendix A. Thus the values for the standard (h s ,z s ) are 
not always well supported. Even so, the effect is substantial and this is an 
important source of systematic error. 

Table 2 is similar to Table 1 in the qualitative sense. The effect of the DV 
profile is greater and at the greater ranges the discrepancies can be quite large. 

4. ERROR ASSESSMENT FOR THREE DIMENSIONAL METHODS 

The ability to locate a sound source position from the transit times (ti, ..., 
t 4 ) needs to be assessed in three dimensions because of the directional 
properties of the array cubes. Our approach is to place the acoustic center at a 
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depth a 2 , and at the center of a right circular cylinder of radius h. The sound 
sources (k in number) are equally spaced on a section of the cylinder at 
depth z. Figure 3 shows a plan view. Azimuth is measured counter 
clockwise with zero at "3 o'clock." 

The ray fitting methodology is used to construct true elevation angles (0i, 
..., 65 ) and true transit times (ti, ..., ts) to each of the k sound sources from the 
four hydrophones and the acoustic center ( 05 , ts). We also need the azimuth 
from the acoustic center (origin of circle) to the sources. These latter are 

<t>j = 2jt(j-l)/k for j = 1, ..., k. (4.1) 

Since the ray fitting methodology is two-dimensional, we must 
characterize the vertical planes connecting each sound source on the cylinder 
to each of the five points in the array cube. The technique for doing this is 
developed in Appendix B. One needs only the horizontal separations and the 
vertical positions of source and receiver. Thus five rays are fit to each sound 
source; one to each of the four phones and one to the acoustic center. Also 
five elevation angles are generated; but we retain only the last, the true 
elevation angle at the acoustic center. 

During operations, the information collected consists of the set of transit 
times (ti, ..., t 4 ) from sound source to receiver array, (see Figure 1). For our 
purposes these four values are regarded as exact. Thus we can use the values 
produced in the ray fitting process. They must be converted to a ray tracing 
direction (azimuth angle, b and elevation angle, 0 ) and a single transit time 
(t ac ) to stop the ray tracing. The currently employed procedure is described in 
[5]. It is convenient to present certain aspects , and this is done in Appendix E. 
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TABLE 1. Comparison of Two Dimensional Ray Tracing. Sensor at 1300 FT. 
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TABLE 2. Comparison of Two Dimensional Ray Tracing. Sensor at 1350 FT. 
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Finally the values (0,t ac ) together with the array depth and DV information 
are fed into the ISOSPEED ray tracing program. The output consists of the 
(h,z) values. Upon combining these with the azimuth 0 C , the estimated 
position of the sound source can be computed. 

A number of error calculations can be made. First the error in estimating 
transit time 

tim er — t a c — ts (4.2) 



affects the rule for stopping the ray tracing algorithm. An error of 10' 4 seconds 
translates to about half a foot in horizontal range, h, (i.e., u = 4800 ft. /sec.). 
Next, the horizontal error is measured directly 

h cr = h c -lM (4.3) 



The values tim er and h er should be strongly correlated. 

In a similar fashion, the error in the elevation angle is correlated with the 
error in the vertical component 



Oer = 0 C ~ 0o 



Z C r — Z c Z] 



(4.4) 



but, because of ray bending, the relationship is non-linear. The value z er is 
measured directly in feet. The value 0 er is more difficult to interpret. It is a 
function of the water layers involved. See [7]. 

Finally we have the error in azimuth 

or 0c 0o- (4.5) 

These values are in radians and, when multiplied by hi, measure how far off 
the mark (see Figure 3) in feet, along the cylinder perimeter. 
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It is a rather incipient fact that our errors are periodic functions of the true 
azimuth. This is illustrated in Figure 4. We suspect that this is at the root 



some earlier attempts to treat possible causes of systematic error, [6,7]. 




DOWN 

RANGE 



Figure 3. Cylinder Cross Section Illustrating Sound Source 
and Receiver Positions 

Many tabulations of these errors have been made. Four sets, each with 15 
points around the cylinder and four radii (hi = 2000, 3000, 40000, 5000 ft.). 
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Seconds Radians Feet 



Figure 4. Periodic Nature of Errors as a Function of Azimuth 

Case: DV of 3/22/89 and Array No. 1. Ref. Table 3. 






have been selected for display and they appear in Table 3. The basis of 
selection was to choose two of the more interesting water columns (10/27/88 
and 3/22/89, see Appendix A) each matched with two of the more severely 
tilted arrays: Number 1 and 56, see Table B-l. The tilts for array number 1 
have opposite signs while those for array 56 are of the same sign. 

Let us make a sample calculation in order to aid in the interpretation of 
these errors. The total error offset is essentially 



Thus for the case of Array No. 1 on 10/27/88 at hi = 5000 ft. and azimuth 
0.418879 radians (i.e.,24 deg. north of east) we have 



in the azimuth, which (measured along the arc) contributes 2.3525 ft. of error 
for every 1000 ft. of horizontal range. A scan of Table 3 shows generally that 
this condition persists, although the periodic nature of the errors can make 
the effect small in some directions. 

It seems that the most severe errors are associated with the arrays having 
the larger tilts. We suspect that the causes lie in the use of the array center as 
the acoustic center and the assumption of constant sound speed for the entire 
array. The z-phone is about 30 ft higher than the others. Also, the tilt 
correction method is but approximate (see Appendix B.) 

We are disinclined to attach great importance to the other source, the 
isospeed ray tracing. The former was considered in Section 3 and cannot 




i^(2. 45) 2 + (.28) 2 + (5-2.3525) 2 =12.02 ft. The dominant portion of the error is 
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Table 3. Error Structure - Current Methodology 

10/27/88 Array #1 Acoustic Center Depth = 1308.76 Sound Source Depth = 250 

hi =2000 0 = 0.483579 t 5 =0.4639118 hi = 3000 0 = 0.334327 t 5 = 0.6521831 

<J) (J) er 0 er tim e r h e r z eT ^er @er tim er h e r z e r 

0000000 .0017878 .0009728 .0000966 -.47 -1.87 .0018177 .0007312 .0000908 -.41 -2.44 
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Table 3. ( Continued ) 
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Table 3. ( Continued ) 
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Table 4. Error Structure - Proposed Methodology 

10/27/88 Array #1 Acoustic Center Depth = 1308.76 Sound Source Depth = 250 

hi = 2000 9 = 0.489930 t 4 = 0.4654736 hi = 3000 0 = 0.339157 U= 0.6532862 
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10/27/88 Array #56 Acoustic Center Depth = 1218.84 Sound Source Depth = 250 

hi = 2000 6 = 0.453929 U = 0.4569741 hi = 3000 6 = 0.312071 t 4 = 0.6473052 
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3/22/89 Array #1 Acoustic Center Depth = 1308.76 Sound Source Depth = 250 

hi = 2000 0 = 0.711452 t 4 = 0.5438750 hi =3000 0 = 0.337225 tj= 0.6551362 



1 oo Os wo OO Os i — i OO oo VO On 
^N^OO(N(N(N(NrHqOO^(N(N 



a 

X 



Os VO i— 

ooooo>— I^OOOOOOO^ 



o 

ID 

<N 

II 

X 



o 

<N 

ID 

CD 

r-H 

ID 

o 



3 



OO SO 
^ (SO 



r^sor^r^oowo^r^coi— <osoo 
^nS-STOCSHOOCSnS- 



D l! 



£ O VO — 

n ^ o o 



^ w <-i 

o o — < 



On cn cs ^ »n On 
o o o o o o o 



0) 

Q 



CD 



oo 






pH 


oo 




CO 




Os 


Os 


r- 


wo 


sc 


VO 




a> 

u 


£ 




N - 


CN 


o 


oo 


Os 


CO 


CO 


wo 


Os 




r- 




VO 


wo OS 


oo 


CO 


N- 


oo 


VO 


wo 


VO 


N- 


CO 


wo 


i—i 


N” 


oo 


VO 




»H 


rt< 




N" 


o 


r- 


n- 




CN 


CN 


N" 


Os 


Os 




vO 


vO 


wo OS 


VO 


VO 


CO 


1— 1 


os 


CN 


vO CO 


Os 


1—1 


CO 


C" 


Os 


vO 


r- 


3 


o 




oo 


r- 


CN 


co 


»— * 


CO 


N* 


vo 


Os 


r- 


CO 


VO 


i—i 


oo 


VO 


oo 


wo 


1—1 


CO 


v£> 


Os 


Os 


oo 


wo 


CO 


o 


CN 


VO 


OO 


Os 


O 


o 


o 


oo 


wo 


i—i 


CO 




Os 


Os 






<n 


i—i 


o 


N - 


r- 


Os 


o 


o 


o 


o 


o 


O 


O 


o 


o 


o 


o 


o 


o 


o 


O 


CD 


<N 


CD 


o 


o 


O 


o 


o 


O 


o 


o 


o 


o 


o 


o 


o 


o 


O 


o 


o 


o 


o 


o 


O 


O 


o 


o 


o 


o 


o 


o 


O 


o 




o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o o 


o 


o 


o 


o 


o 


O 


O 


o 


o 


o 


o 


o 


o 


O 


o 


T3 


O 




o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


O 


o 


o 


o 


o 


o 


o 


o 


o 


o 


£ 


II 




o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 



o 

CD 



CD 



CO 


wo 


Os 


o 


VO 


oo 


CO 


VO 


N- 


VO 


CN 


VO 


wo 


co 


N- 




o 




wo 


CN 


VO 


CO 


1—1 


oo 


wo 


CO 


N- 


Os 


VO 


CO 


N- 


o 


VO 


vO 


Os 


oo 


CO 


oo 


VO 


Os 


r— t 


CN 


1—1 


Os 


OO 


N- 


N* 


00 




o 




oo 


VO 


1—4 


o 


r-* 


N- 


CO 


CN 


N- 


VO 


CN 


CO 


CN 


co 


wo 




CO 


Os 


oo 


Os 


VO 


1—4 


i-H 


o 


wo 


vO 


»o 


o 


Os 


wo 




o 


(m 


CO 


1—4 


wo 


N“ 


oo 


o 


pH 


VO 


o 


1-H 


CN 


CO 


CN 


VO 


CO 


CN 


CO 


CO 


CO 


CN 


1—4 


o 


1—1 


CN 


CN 


CN 


CN 


CN 


o 


o 




ID 


o 


1-H 


CN 


CN 


CN 


1—4 


1-H 


o 


o 


1-H 


1-H 


1-H 


1—4 


i-H 


o 


o 


O 


o 


o 


o 


o 


o 


o 


o 


o 


o 


O 


O 


o 


o 


o 






o- 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 




II 




o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


vO 






o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


CN 






o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 



00 

O 

CO 



hrl-OsCNWOi— 'WOSO tJ- 
^ ^ O O O H rH H 



o ^ »o 

CO CO 1—4 ^h 1—4 



J(N00(N^Of0^r)O>OO^On^ 
^ ^ O O p O ^ ^ ^ CD O O O »-« ^ ,-h 



Oh O 

0> CN 5nwocow“>cooor^oocNCNCNcot^^coos 
Q 00 tsfcocNo^cNcocococNi-<oocNcoco 

^ iii i i i i 



a> 

•4-# 

o 

a> 

U 

o 



ID 

oo 



- X* 

4-» 

CO 
0 



<u 

x 



0\O^S”COO 

1— • o o o o —• 



»— 'OvsOcO^CNsoON 
i— 'OOOOOOO 



O 

CD 



N- 


wo 




VO 


wo 


CO 


CO 


CO 


wo 


CO 


Os 


CO 


CO 


Os 


Os 


o 

o 


VO 




r— 4 


CN 


CO 


oo 


1— < 


N- 




1— < 


CO 




o 


wo 








VO 


wo 


VO 


N" 


VO 


CO 


N’ 


Os 


CO Tfr 


o 


r- 


wo 


r- 


1-H 


< 






oo 


oo 


o 


oo 


Os 




VO 


1—1 


o 


Os 


VO 




r- 


OO N" 


oo 


N- 


1—4 


CO 


CO 


CO 


Os 


1-H 


CO 


r- 


CN 


1—4 


r- 


N - 


Os 


^4 


CD 




r- 


so 


co 


CN 


o 


CO 


wo 


o 


wo 


Os 


r- 


oo 


1-H 


CN 




VO 


N - 


1—4 


CN 


wo 




r- 


r- 


wo 


CN 


o 


CO 


wo 




r- 




LD 


qj 


oo 


wo 


1—4 


CO 


r- 


Os 


Os 


oo 


wo 


CN 


o 


1-H 


wo 


oo 


Os 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 




CN 


CD 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o o 


o 


o 


o 


o 


o 








o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 




o 




o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 




II 




o 


o 


o 


o 


o 


o 


o 


o 


o 


o o 


o 


o 


o 


o 



CD 



=tt= 



wo 


o 


CN 


Os 


Os 


Os 


CO 


CO 


oo 


CO 


o 


oo 


CO 


wo 


r- 






o 


1—4 


o 


oo 


N- 


Os 




VO 


oo 








o 


VO 


CO 


oo 


o 


o 


VO 


r- 


VO 


N- 


VO 


1-H 


VO 


oo 


1-H 


r- 


l> 


CN 




o 




oo 


CO 


oo 


N- 


os 




wo 


o 


o 


C-* 


oo 


o 


CN 


CO N - 


Os 


CN 


CO 


o 


wo 


1—4 


CO 


co 


wo 00 


CN 


r- 


1-H 


r- 


1-H 


Jh 


o 


»h 


VO 


VO 


o 


Os 


CN 


CN 


1-H 


oo 




SO 


r- 


oo 


W-> 


r- 




CO 


VO 


r- 


r~ 


wo 


CO 


o 


CN 


N* 


wo 


VO 


wo 


N" 


1-H 


1—4 


Jh 


rf 


QJ 


r-H 


CN 


CO 


CN 


CN 


1—4 


o 


o 


1—4 


i-H 


f— * 


1-H 


1—4 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


< 




-o- 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


^4 


II 




o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 




l 




o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o o 


o 


o 


o 


o 


o 




X 




o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


O p 



o 


o 


o 


o 


o 


1-H 


1-H 


1-H 


1-H 


1-H 


1-H 


o 


o 


o 


o 




o 


o 


o 


o 


o 




1—4 


1—4 




1-H 


1—4 


o 


o 


o o 


o 


Os 


oo 


r- 


VO 


wo 




CO 


CO 


Of 


wo 


VO 


r- 


oo 


Os 




o 


Os 


oo 


r- 


SO 


w^ 




CO 


CO 




wo 


VO 


r- 


oo 


Os 


o 


r- 


wo 


CO 


1—4 


Os 


r- 


wo 


wo 


r- 


Os 


i-H 


CO 


wo 


r- 


o 


r- 


wo 


CO 


1—4 


Os 


r- 


wo 


wo 


r- 


Os 


1—4 


CO 


wo 


r- 


o 


oo 


r~ 


VO 


wo 


CO 


CN 


1—4 


1-H 


CN 


CO 


wo 


so 


r- 


oo 




oo 


r- 


VO 


wo 


CO 


CN 


1-H 


1—4 


CN 


CO 


wo 


VO 


r- 


oo 


o 


oo 


r- 


VO 


wo 


O}- 


CO 


CN 


CN 


co 


N* 


wo 


SO 


r- 


oo 


CD 


^ o 


oo 


r- 


VO 


wo 




CO 


CN 


CN 


CO 




wo 


VO 




oo 


o 


1-H 


CO 


wo 


r- 


Os 


1—4 


CO 


CO 


1—4 


Os 


r-* 


wo 


CO 


1—4 


o 


1—4 


CO 


wo 


r- 


Os 


1—4 


CO 


co 


1-H 


Os 


r- 


wo 


CO 




o 


p 


oo 


CN 


VO 


p 


wo 


Os 


p 


wo 


p 


p 


CN 


oo 


p 




o 




oo 


CN 


vq 


p 


wo 


Os 


ON 


wo 


p 


vq 


CN 


oo ^ 








f— < 




CN 


CN 


CN 


CN 


CN 


CN 




X 


1 


1 










1— < 


f— < 


CN 


CN 


CN 


CN 


CN 


CN 


1—4* 




1 


1 



i i i i i i i i i i 



32 



3/22/89 Array #56 Acoustic Center Depth = 1218.84 Sound Source Depth = 250 

hi = 2000 0 = 0.452549 t 4 = 0.4583009 hi = 3000 8 = 0.309991 t 4 = 0.6491808 
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support errors of this magnitude. (The calculations in this section apply the 
polynomial extrapolation prior to the development of isospeed layers. In this 
way the DV table extension is not an issue in the present comparisons.) 

Accordingly, the author has proposed and tested a method that treats the 
two main issues. The acoustic center is moved to the C-phone. (Visualize 
Figure 3 with the C-phone on the axis of the cylinder.) In this way the transit 
time that stops the ray tracing is a directly measured quantity. Also the 
sound speed in the layer containing the (X, Y, C) phones (Z-phone omitted) is 
assumed constant in order to determine (<|>p, 0 p ), the azimuth and elevation 
angles of the proposed method. In this way the change in sound speed at the 
Z-phone is taken out of the computation. Details appear in Appendix F. 

Computations using this technique appear in Table 4. The improvement 
is dramatic. Let us repeat the earlier exemplary computations 10/22/88, Array 
No. 1, hi = 5000 and <{) = 0.418879. This time the inputs come from Table 4: 
specifically V (.22) 2 + (0.5) 2 + (5-0.04306) 2 = 0.31 ft., a value dramatically smaller 
than 12 ft. A scan of Table 4 and comparison with Table 3 shows that these 
improvements are quite consistent. 

These computations also have the advantage that the five-foot layer 
isogradient ray tracing was used together with the exact tilt corrections, eq. 
(B.5.) 

5. CONCLUSIONS 

Generally there are a number of sources of systematic error. In some the 
effects are small, e.g., isogradient vs. isospeed raytracing, 5 foot versus 25 foot 
layer thicknesses. None the less they are systematic and certainly no longer 
necessary. The presence of systematic errors tends to build up idiosyncracies 
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that frustrate the use of standard statistical methodology when 
troubleshooting other aspects of the data. 

The discovery of the periodic nature of the errors was a surprise. Its 
presence adds another dimension to the interpretation of the data. No doubt 
it is the source of some deception. 

The major sources of error related to ray tracing (see pages 5, 6) are 
believed to be associated with 

(i) the conversion of transit times to the direction of the sound 
source 

(ii) the constant speed extrapolation of the DV profile below 1300 ft. 

(iii) the use of approximations for tilt corrections. 

Of course the effects of these errors are directional because of their periodic 
nature. We have documented errors of 12 ft. or so due to (i). We can 
speculate 10 or more feet because of (ii). Moreover the combined effects could 
be additive. The effect of (iii) increases as the tilts increase. 

It is further noted that the above errors apply to single arrays. At this 
point we have no comment about how they may combine to produce 
mismatches in the array overlap regions. It would be more comfortable to 
treat this issue by introducing the changes and then collecting more data. 

6. ADDENDUM 

It came to the author's attention during the final editing phase of this 
report, that the approximate tilt corrections, eq. (B.4), are not the ones 
currently employed by the software. Rather, the system is using 



X(l) 




c 2 0 -S 2 


x 0 (l) 


X(2) 


► = < 


0 c-i —Si • 


• x 0 (2) 


X(3) 




S 2 Si CiC 2 


X 0 (3) 
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Accordingly the author made some specialized computations for 
purposes of indicating the effect. Referring to Tables 3 and 4, we have chosen 
to treat the exemplary case: water column of 10/27/88, array number 1, h] = 
5000 ft, zi = 250, and <}> = 0.418879 radians. The table below contains the four 

errors and total offset for each of four algorithms: 

(i) the isospeed method (from Table 3) 

(ii) the isospeed method using (6.1) vice (B.4) 

(iii) the isospeed method using (B.5) vice (B.4) 

(iv) the proposed method (from Table 4) 

TABLE 5. EFFECT OF TILT CORRECTION METHOD 

10/27/88, Array #1, h] = 5000 ft, Z\ = 250 ft, ((> = 0.418879 rad., a 2 = 1308.76 ft. 





4*er 


0 e r 


her 


Z e r 


d 


(i) 


.0023525 


.0003869 


.28 


-2.45 


12.02 


(ii) 


.0021174 


-.0001567 


1.47 


2.95 


11.09 


(iii) 


-.0000161 


-.0000077 


1.26 


2.19 


2.53 


(iv) 


.0000166 


.0000430 


-.05 


-.22 


0.31 



The results of Table 5 suggest the following: there is little distinction 

between the use of (B.4) and (6.1) for the tilt correcting (compare (i) and (ii)). 
However case (iii) suggests that there is much to be gained by use of the exact 
tilt correct (B.5) even if isospeed ray tracing is used with 25 ft layer thickness. 
Finally case (iv) suggests that considerable gains are available if, in addition, 
we use the proposed methodology. 

All of these systematic errors are mathematical, i.e., due to choice of 
algorithms. There is no longer any reason not to use the best. 
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APPENDIX A 



The twelve depth velocity profiles used in this study are recorded here in 
graphical form. They include the profiles used in [7]. Generally the measured 
values stop at a depth of about 1300 feet and it is necessary to extrapolate, in 
several instances to depths greater than 1350 feet (see Table B-2). 

The insets of the graphs illustrate two different extrapolation schemes. 
The constant value extrapolation is the one currently in use, [5]. But there is a 
slight deception. Current methodology utilizes isospeed profiles, see Figure 2, 
and 25 foot layer thicknesses. The constant value extrapolation is the value of 
the deepest 25 foot layer appearing in an isospeed profile. 

The other extrapolation is based upon fitting a second order polynomial 
to the deepest hundred feet of the original profile (five foot layer thicknesses). 
There are two steps in this process. First is fitting the curve by least squares. 
Because of the equally spaced depth increments, the fitting takes an especially 
simple form. Using the equation, 

v = a + b(n-u) + c(u-u)^ (A.l) 

where x> is velocity, u is layer depth, and u = average depth, the normal 
equations take the form 

Xv = Xa - Xb(u-u)- cX(u - i7) 2 
Xv(u-u) = Xa(u -17)- Xb(u-u) 2 -cX(u-u)^ 

Xv(u-u ) 2 = Xa(u-u) 2 ' - Xb(u-uf > -cX(u-u) 4 
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Using the notation S v = £v, Sx> u = £v (u- u ), S^uu = £v (u- u ) 2 . S U 2 = £ (u- 
u ) 2 , S u 4 = Z(u- u ) 4 and recognizing E(u-u)= Z{u-uf > - 0 because of the 
uniform spacing, the above equations assume the reduced form 

Sv= na + cS u2 

= bs„ 2 (A. 2 ) 

S vuu ~ aS u 2 cS u 4 

and we solve for the coefficient of the quadratic term 

c - (nS vuu — S V S U 2) / 4 — S u 2 j (A. 3) 

(This done, values for a and b come easily from the first two equations.) 

We note in passing that if n = 2K+1, an odd number, and A is the spacing 
between consecutive values of (u), then u is an integer and 

2 

S„2=yK(K + l)(2K + l), 

s«4 = y[ k ( k+1 )] 2 ' <A ' 4) 

nS u 4 - S^ 2 = A 4 K 2 (K + 1) 2 (2 K + 1 )(7 - 4K) / 18 

Such formulas expedite the computation of c in (A.3). 

The second step deals with the issue of how to use this quadratic for 
extrapolation purposes. Generally, the direct use of the equation (A.l) would 
lead to a visual discontinuity between the last measured value and the first 
extrapolated one. We have chosen not to do this. Instead we recognize that 
A c represents the first difference of the gradient sequence {x>j (i)}, see Figure 2. 
So to preserve continuity, the extrapolation is enabled by using successive 
updates 
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L»i(l) = l>[ (/ — 1) + Zi -C 

u(i + 1) = v(i) + Avi(i) 
v 0 (i) = [u(i + l)u(i) - u(i)v(i + 1)] / A 

where \)(i) is the estimated sound speed at u(i); \)o(i) and \>i(i) are the intercept 
and slope values of the straight line fit for the i th layer. It is this continuity 
preserving step that explains the occasional appearance of "crooked" 
extrapolations in the insets. 

The use of isogradient ray tracing with 25 foot layer thicknesses was used 
in Section 3. If extrapolation of the water column was required, then a 
different second order polynomial method was used, and after the conversion 
to 25 foot layers. Basically the quantity A c was estimated by averaging the 
difference of the last five values of \)i (last 100 feet), and then proceeding in 
the same way as stated above. No graphs showing the effect of this have been 
prepared. But the choice does have an effect upon the d g values computed for 
Table 2. 

We take this opportunity to note that the visual appeal of the quadratic 
extrapolation is rather good in instances 1, 2, 3, 4, 6 and 11. The others are 
easily challenged. This information may influence the reader's interpretation 
of some values appearing in Table 2. 
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Depth Velocity Profile — 5/12/88 
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Depth Velocity Profile — 6/22/88 



41 




4860 4880 4900 4920 4940 

Velocity ( ft/sec ) 

Depth Velocity Profile — 7/21/88 
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Depth Velocity Profile— 8/03/88 
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Depth Velocity Profile — 10/27/88 
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Depth Velocity Profile— 1/13/89 



45 



r 



T 



T 




Depth Velocity Profile— 3/08/89 
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Depth Velocity Profile — 3/22/89 
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Depth Velocity Profile — 4/26/89 
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Depth Velocity Profile— 4/27/89 
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Velocity ( ft/sec ) 

Depth Velocity Profile — 5/10/89 
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Depth ( ft ) 
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Velocity ( ft/sec ) 



Depth Velocity Profile — 6/06/89 
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APPENDIX B 



COORDINATE SYSTEMS 

The developments that follow deal with several coordinate systems and 
we need efficient ways to distinguish among them. First of all is the right 
handed system defined by the X, Y, Z and C hydrophones of the array. All 
incoming transit times must be interpreted in this system and we call it cs(a), 
or the coordinate system of the array. The origin is at the c hydrophone. 

Since this system is generally tilted with respect to a "flat earth surface" 
system there is need to rotate cs(a) into alignment with a common coordinate 
system for all arrays, or a range coordinate which has horizontal directions 
consistent with the earth's surface. Such a resultant system will be called cs(b) 
and, for convenience of terminology, the X arm of cs(a) is rotated to a position 
called east; the Y arm to a position called north, and the Z arm to a position 
called vertical or zenith. The origin is still at the c hydrophone. 

The ray tracing methodology of Ref [5] attempts to locate the sound source 
in relation to a specific point termed the acoustic center. This center is the 
geometric center of the array cube. The resulting coordinate system is a 
translation of cs(a), and will be called cs(ac). In this system the four 
hydrophone positions are specified by D/2 times the vectors 
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respectively, whereas in cs(a) these four vectors would be the three unit basis 
vectors and the origin. Conversion from cs(ac) to cs(a) is a direct translation. 

The conversion of a point located in cs(a) to the cs(b) system requires a 
three-dimensional rotation based upon the tilt angles and the "ZROT" 
horizontal direction correction. It is convenient to describe this in terms of 
the three Euler angles <t>i, <J> 2 , <1>3, or roll, pitch and yaw. Letting 

Si = sine (4>i) and Ci = cosine (<t>i) i = 1, 2, 3 (B.2) 



we define three successive rotations which, when applied sequentially to a 
point in cs(a), will in the end describe it in cs(b). First hold the X arm fixed 
and rotate the Y-Z plane through an angle <J»i; the matrix of this transform is 



Pi 



1 0 0 

0 C] -si 

- 0 si ci _ 



Next hold the (current position of) Y-arm fixed and rotate the X-Z plane 
through an angle of <{> 2 ; this transformation has matrix 



C2 0 -52 



P2 = 



0 1 
_ S2 0 



0 

C2 - 



Finally hold the (current position of) Z-arm fixed and rotate the X-Y plane 
through an angle of (fo; the transform is 



C3 


-S3 


S3 


C3 


_ 0 


0 



0 

0 

1 _ 



The successive applications of these three rotations is a unitary 
transformation, (i.e., its inverse is equal to its transpose). 
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B = p3 p2pl 



(B.3) 



and if a is a three vector in cs(a), then b = Ba is the same vector referenced in 
cs(b). 

The determination of the Euler angles is accomplished as follows. The 
submerged arrays have tilt indicators on the X and Y arms which, 
individually, measure the angles that these arms make with the horizontal. 
An accounting for these tilts must be made when the ray trace azimuth and 
elevation angles are converted to a horizontal based coordinate system. The 
technique currently in use takes the apparent position, Xo, and applies the 
transformation 



'X(l)' 




1 


0 


-sin(XTILT) 


'Xo(l)' 




•X(2) 


> — < 


0 


1 


-sin(YTILT) ► 


• X 0 (2) 


• (BA) 


X(3). 




sin(XTILT) 


sin(YTILT) 


1 


X 0 (3), 





so that the new apparent position, X, is referenced in a plane level with the 
earth. This transformation is an approximation which simply rotates the two 
arms to the horizontal as if they were separate unconnected arms and the 
rotation of one does not affect the rotation of the other. That is, the first two 
rows of the coefficient matrix are not orthogonal. The result is an 
approximation whose success depends upon the smallness of the tilt angles. 

The exact way to accomplish this goal involves the direct replacement of 
the coefficient matrix with the product rotation p2pT 



X(l) 




^2 


-S]S 2 


-C]S 2 ' 




'Xo(l) 


X(2) 


> = < 


0 


C 1 


-S 1 




Xo(2) 


X(3). 




. S 2 


S]C 2 


c,c 2 




X„(3) 
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Upon comparing these two coefficient matrices, one sees there is choice in 
identifying the two tilt angles with these two Euler angles. We have chosen 
to match the first two elements of the third row: 

s 2 = sin(XTILT) s i = sin(YTILT) / c 2 . (B.6) 

The geometric interpretation is as follows. First hold the X arm fixed and 
rotate the plane of the Y-Z arms so that the Y arm is horizontal. This is not a 
vertical projection. The division by c 2 shows that one must rotate through an 
angle greater than YTILT in order to maintain orthogonality of the coordinate 
system when making the Y arm parallel to the earth's surface. This done, we 
next hold the Y arm fixed in its new position and rotate the plane of the X-Z 
arms so that the X arm is horizontal. Since the new Y arm is already 
horizontal this second rotation is a vertical projection through an angle of 
XTILT. 

This latter method is exact. The nature of the original approximation can 
be assessed by comparing the two coefficient matrices, using numerical 
inputs. The effect is not great for most of the tilts present at Nanoose. 

To complete the conversion of cs(a) to cs(b) we must find the third Euler 
angle in terms of ZROT, the rotation of the horizontal plane determined by 
array survey. In [ 5 ], ZROT is defined as the angle from the horizontally 
projected X arm to the range center line (east). The comparison of this 
definition with that of <J>3 (see P3) leads to 

s 3 = -sin (ZROT) (B. 7 ) 

The subsequent application of p3 to P2P1 will bring the point a into east-north 
orientation. 
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Finally the position location system prefers to specify an object's position 
in terms of east, north, and depth (positive) below the sea surface. This 
system is a left-handed one with origin on the sea surface directly above the 
acoustic center of the array. 

Table B-l contains the positions of the Nanoose arrays at the time of their 
most recent survey dates. These are the positions of the acoustic centers in 
the range coordinate system. It is both useful and instructive to use our 
coordinate system superstructure and locate the hydrophones of each array in 
the range coordinate system. 

Let a be the 3- vector locating the acoustic center of an array in the range 
coordinate system, see Table B-l. Let a be the location of one of the phones in 
cs(a) and let f, 



m 



f = (D/2) 



(B.8) 






be the position of the acoustic center in cs(a). In the system cs(b) these vectors 
are Ba and Bf, and the location of the phone relative to the acoustic center is 
B(a-f). The position p of the phone in the range coordinate system is 

p = a+B(a-f) (B.9) 



Since the C-phone is the origin in cs(a), its range coordinate position is a - Bf. 
the result of this computation is in Table B-2. 
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TABLE B-l. NANOOSE ARRAY COORDINATES AND ORIENTATION 

ANGLES 

(Distances in feet; angles in radians) 



Survey Date 


X 


Y 


Z 


XTILT 


YTILT 


ZROT 


ARO 


6/20/85 


12188.01 


-131.52 


-1295.33 


0.002909 


0.014835 


-0.208183 


AR1 


6/20/85 


19463.16 


-174.99 


-1308.76 


0.061523 


-0.036070 


1.362579 


AR2 


7/12/85 


26991.39 


-109.83 


-1323.25 


0.000145 


0.005236 


2.670336 


AR3 


1/7/88 


34505.10 


-80.76 


-1323.32 


0.027925 


-0.011345 


2.928139 


AR4 


10/24/88 


42005.19 


-55.17 


-1318.28 


0.001164 


-0.040288 


-2.315877 


AR 5 


6/20/85 


49497.00 


-25.23 


-1315.58 


-0.000291 


-0.004072 


1.668535 


AR6 


6/20/85 


56972.28 


-21.21 


-1308.50 


0.013817 


0.041161 


-0.703420 


AR7 


7/30/85 


64680.66 


15.33 


-1353.39 


0.034907 


0.022835 


-0.574144 


AR8 


11/16/88 


71969.73 


-29.28 


-1300.89 


-0.005963 


-0.012217 


-1.577341 


AR9 


5/7/84 


3.00 


3.00 


1.00 


0.000000 


0.000000 


0.000000 


AR 10 


3/12/84 


47100.00 


-3600.00 


-1300.00 


0.000000 


0.000000 


0.000000 


AR 11 


7/18/85 


23173.89 


-6488.40 


-1312.09 


-0.004654 


0.000436 


2.784376 


AR 12 


6/20/85 


30731.25 


-6553.05 


-1312.90 


0.002036 


0.001745 


-3.042179 


AR 13 


6/20/85 


38213.61 


-6640.77 


-1323.05 


0.000291 


0.006254 


1.373522 


AR 14 


6/20/85 


45647.07 


-6513.18 


-1324.78 


0.001309 


0.002327 


-2.348044 


AR 15 


6/19/85 


53249.43 


-6354.60 


-1316.66 


0.003345 


0.004509 


0.581544 


AR 16 


9/13/85 


60859.74 


-6356.07 


-1313.42 


0.014835 


0.036943 


2.303276 


AR 17 


6/16/87 


68217.93 


-6524.10 


-1313.43 


0.008290 


0.034761 


2.158449 


AR 54 


2/2/88 


38029.95 


5401.98 


-1212.69 


0.007709 


-0.003782 


-1.056919 


AR 55 


6/20/85 


45645.75 


6369.66 


-1188.12 


0.027634 


0.039415 


-0.728553 


AR 56 


7/30/85 


53180.13 


6417.96 


-1218.84 


0.037525 


0.048142 


-1.392651 


AR57 


7/30/85 


60745.71 


6419.40 


-1088.24 


0.006981 


0.001891 


-3.108606 


AR 23 


6/20/85 


41605.14 


-12150.18 


-1268.23 


-.002182 


0.003200 


-1.845214 


AR 24 


4/17/89 


49572.00 


-12966.00 


-1300.00 


-0.007272 


0.055269 


-1.343904 


AR 25 


10/24/88 


56993.79 


-12999.33 


-1205.48 


0.000291 


-0.002182 


-0.593726 


AR 26 


8/8/88 


64442.94 


-12971.04 


-1255.35 


-0.014835 


-0.012654 


3.134192 


AR 27 


7/15/80 


22119.60 


-15908.70 


83.00 


0.000000 


0.000000 


0.000000 


AR 28 


5/4/83 


45000.00 


1500.00 


-1350.00 


0.000000 


0.000000 


0.000000 


AR 29 


2/2/79 


0.00 


0.00 


0.00 


0.000000 


0.000000 


0.000000 
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AR 

0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

54 

55 

56 

57 

23 

24 

25 

26 

27 

28 

29 



TABLE B-2. LOCATIONS OF THE C-HYDROPHONES 



Date 


xc 


YC 


zc 


6/20/85 


12170.190 


-143.316 


-1310.105 


6/20/85 


19473.791 


-192.189 


-1325.075 


7/12/85 


27011.563 


-103.354 


-1338.287 


1/7/88 


34522.511 


-69.103 


-1338.681 


10/24/88 


42004.318 


-33.398 


-1332.413 


6/20/85 


49513.397 


-38.633 


-1330.630 


6/20/85 


56950.933 


-23.551 


-1323.123 


7/30/85 


64659.408 


10.557 


-1367.552 


11/16/88 


71954.918 


-13.999 


-1315.793 


5/7/84 


-12.000 


-12.000 


-14.000 


3/12/84 


47085.000 


-3615.000 


-1315.000 


7/18/85 


23193.258 


-6479.598 


-1327.004 


6/20/85 


30744.657 


-6536.662 


-1327.956 


6/20/85 


38225.375 


-6658.513 


-1337.943 


6/20/85 


45646.877 


-6492.002 


-1339.829 


6/19/85 


53245.085 


-6375.441 


-1331.552 


9/13/85 


60880.699 


-6357.757 


-1328.681 


6/16/87 


68238.605 


-6528.792 


-1328.448 


2/2/88 


38009.399 


5407.725 


-1227.510 


6/20/85 


45624.165 


6367.887 


-1202.470 


7/30/85 


53162.159 


6429.360 


-1233.742 


7/30/85 


60760.102 


6434.858 


-1103.370 


3/20/85 


41594.799 


-12131.725 


-1283.312 


4/17/89 


49554.120 


-12955.612 


-1315.728 


10/24/88 


56972.961 


-13003.338 


-1220.483 


8/8/88 


64458.271 


-12955.966 


-1269.935 


7/15/80 


22108.374 


-15890.700 


68.000 


5/4/83 


44985.000 


1485.000 


-1365.000 


2/2/79 


-15.000 


-15.000 


-15.000 
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APPENDIX C 



2. RAY FITTING AND RAY TRACING 

A "ping" sound source produces a wave front that travels through the 
water and is detectable by the receiving transducers. A ray is the path 
generated by the normal to the wave front. Ray tracing is the activity of 
following a ray from a receiver for a fixed amount of time in order to locate 
the sound source. Ray fitting is the activity of recreating the ray path from the 
positions of the source and the receiving sensor. The latter is needed to study 
the error performance of the former. 

The speed of sound in water is assumed to vary with depth, but remain 
homogeneous in the horizontal. Thus a single depth-velocity profile is valid 
throughout the field. 

We begin with some preliminaries followed by the development of a ray 
fitting algorithm, which may be viewed as an inverse to ray tracing. It is 
certainly more difficult. All paths are direct paths. That is, no provision is 
made for reflections or refractions that produce non monotone rays 
connecting source and receiver. Also our interest lies in the greater ranges 
and no adjustments are presented for sources that are directly above the 
receiver. 

The notation is chosen to be consistent with that used in the Fortran 
source codes. We are given a speed-depth profile in the form of pairs lm, and 
velj. 

lmj = depth of ith water layer, positive down, 
velj = speed of sound at the depth lm,. 
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Digression. The water speed processing system at NUWES produces 

average velocity values for a large number of equally spaced points. These 

averages are intended to serve as constant values for the entire layer, eq. (3.1). 

Thus the information is provided in the form, for i = 1, ..., m 

li = depth of water layer boundary 
velj = velocity constant for layer (lj, l, + i) 

The values used for lmj above are the layer midpoints: 
lm; = .5 x (lj + lj+i) for i=l, ..., m-1 

The algorithms developed below do not require layers of equal thickness. 
Thus they can accommodate the user who wants to use thin layers at depths 
of rapid change and thick layers at depths of slow changes. All computations 
should be in double precision arithmetic. 

Occasionally as application leads to a receiver whose depth is greater than 
the largest {lmj}. For these cases we have been using an extrapolation of the 
sound speed profile that adjoins a sufficient number of layers, all of the same 
thickness as the deepest of the original layers. The corresponding velocities 
are extrapolated using a second order polynomial calculated from a fixed 
second difference whose value is the average of the four deepest second 
differences. (Use of the coefficient of the quadratic term in a least square fit 
has also been employed as an option.) 

Since ray fitting is a rather delicate operation we use the isogradient 
technique. I.e., straight lines are fitted to the speed for each layer; the constant 
gradient of slope is computed; the profile itself is then a continuous function 
of connected straight line segments. See Figure 2. I.e., if 

lmj < z < lmi+i 
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then 



vel(z) = vo(i) + vi(i)z 

where 

vo(i) = lntj+i • velj-lmi • veli+i 
vi(i) = (vel i+ i-veli)/dzi 
dzj = lm i+j-lmj 

Typically the values of {vo(i)} are positive and large compared to the 

{ Vi (i)} which are small and of either sign or zero. Snell's law comes into play 

and the ray invariant will be denoted rv. The ray path in each layer will be a 

circle arc (vi(i)*0) or a straight line (v 1 (i)=0). 

Now we are positioned to describe our ray fitting algorithm. It is two 

dimensional (horizontal, vertical) and given the endpoints 

(ai,a2) receiver 
(PI/P2) sound source 

the goal is to compute 

00 entrance angle at the receiver 

01 exit angle at the sound source 

t transit time of sound from source to receiver 

There is no loss in placing the origin on the water surface directly over the 
receiver. Thus, a\ = 0 and p2 < a2 (depth of receiver). 

The algorithm is an iterative one that operates as follows. Initialize the 
entrance angle 0o- Use this angle to ray trace upward through the layers until 
the vertical value of p2 is achieved. Compute the current horizontal value h. 
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and compare it with pi. If this is within a preassigned small number, e, stop 



Figure C-l. Details of single layer processing 



v(z) = vo + viz (linear profile) Zi<z<zo 
Given: 6 0/ z 0 , z h v 0 , vi, and 



the ray invariant rv = ~ 

Case V] = 0 
dz = Z] — zq 
dw = dz/sin(0o) 
h] = ho + dw • cos(0o) 
dt = dw/vo 
cos(0i) = rv • v(zi) 

Case vi * 0 

q2 = -vo/vi 

s = sin(0o) 

c = cos(0 o ) 

qi = ho + (q 2 - z 0 )s/c 

r = signum(q 2 )(q2 - zo)/c 

h = q, - signum(q 2 ) r • sin(0) 



cos(0o) 




t>, ^ cos(0) Je 0 



1 . ( l + sinf#)^ A 
— In — ^ / 
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co;(0i) = rv • v(zj) 



and compute the transit time. Otherwise, adjust the entrance angle 0 q 



h 



a 2 -p 2 



according to the ratio of the current rise over run. 

Pi 

a2 _p 2 ' and repeat the algorithm using the new initialization. 



-, to the desired one. 



ALGORITHM RAYFIT 
Initialize 

(i) Determine the layers that contain the source and receiver. 
Choose j, n so that 



lmj < p 2 < lm j+1 



lm n < a 2 < lm n+ i 

(ii) Make thickness corrections for the extreme layers 



dZj = lm j+1 -p 2 



dz n = a2~lm n 

(iii) Compute the sound speed at depths a2, p 2 

va 2 = v 0 (n) + vi(n) • a 2 



Vp2 = V 0 (j) + V](j) • p2 

(iv) Unless a "previous" value for 0o is available, fit a straight line 
through the depth- velocity profile in the range (p 2 , a 2 ) and use 
00 corresponding to the circle arc (or line) of that approximate 
profile. I.e., 

If va 2 = vp2 then 0o = tan -1 ((a2 = p2)/pi), else 

va 2 • p 2 -vp2 ■ a 2 

q2 “ (a 2 -p2> 



qi = .5 • pi + .5 • (p 2 -a2)(p2+a2-2q 2 )/pi 
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00 = tarr 1 {qi/(q2— a 2 )} 



endif 



Set initial values for iteration 

A. i = n, 

rv = c/ va 2 . 

Main raytracing code 

B. If vi(i) = 0, then 



s = sin(0 o ), 
h = 0, 



c - Vl-s 2 
z = a 2 



dw = dzj/s 
h = h + c • dvv 



else 



q2 = -vo(i)/vi(i) 
qi = h + (q 2 - z) • s/c 

r = a/ (h-qi) 2 + (z-q 2 ) 2 
c = rv • vel(i) 

s = V 1— c 2 

h = qi - signum(q 2 ) • r • s 



z = lnij 
i = i-1 
goto B 

TEST: 0i = cos -1 (rv • vp 2 ) 

If vi(j) * 0 h=qi-signum(q 2 ) • r • sin(0i) 
If I h-p 2 1 < e , goto FINI 

Re-estimate 0o 

0 O = tan _1 {tan(0o) • h/p^ 
goto A 
FINI: 

ang(j) = 0; 
ang n+ i=0o 

angi = cos _1 (rv • vel(i)) for i = j+1, n 



endif 



If i = j. 



goto TEST 
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Compute transit time 
t = 0 

Do TC i = j, n 
If vi(i) = 0 

t = t + dZj/(vo(i) • sin(angi)) 

else 

_ 1 fcos(ang i+ i)(l+sin(angj)) | 

1 “ vi(i) ln j(l+sin(ang i+1 ))cos(angi)J 
endif 

TC continue 

Remove extreme layer thickness corrections 
dzj = lmj + i - lmj 
dz n = lm n+ i - lm n 

end 

This algorithm has been quite useful to us. Using e = 10~ 6 we typically 
have 8 to 10 iterations through A. 

Raytracing algorithms are less sensitive, but since a good one is readily 
available by merely modifying the above, let us do that. The process is an 
inverse one in that the goal is to compute pi, p 2 given 0o and to, where to is 
the transit time. 

ISOGRAD 

Initialize by locating the layer containing the receiver, establishing the ray 
invariant, etc. 

Choose n so that lm n < a 2 <lm n+ i 

i = n, h = 0, z = a 2 , t = 0 

s = sin(0o), c = cos(0o) 

va 2 = vo(i) + vi(i) • a 2 

rv = c/va 2 

dz n = a 2 - lm n 
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A A: If Vj(i) = 0, then 



dw = dzj/s 
dt = dw/ vo(i) 
h = h + c • dw 

else 

q2 = -vo(i)/v 1 (i) 
qi = h + (q 2 -z) • s/c 

r = V(h-qi) 2 + (z-q 2 ) 2 
cp = rv/vi(i) 



sp = Vl-cp 2 




h = qi-signum(q 2 ) sp 



endif 

t = t + dt 

If t > to goto FINAL 
z = lrrii 

s = sp c=cp i = i-1 goto AA 

FINAL 

dt = t 0 + dt - t 
If vi(i) - 0 
dzj= v 0 (i) • dt 
dw = dzj/sp 
pi = h + dw • cs 
P 2 = z- dZj 

else 

x = exp{dt • V](i)}(l+s)/c 
cp =2 • x/(l+x 2 ) 

sp = V 1-cp 2 
P! = qi + r • sp 
p 2 = q 2 + r • cp 
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endif 

Restore layer integrity 



dzj = lmj+i - lmj 
dz n = lm n+] - lm n 

Compute exit angle 

0i = cos _1 (cp) 
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APPENDIX D. EFFECT OF SOUND SPEED OSCILLATION WITHIN A 

LAYER OF WATER 

Suppose a water layer of thickness A z has an average sound speed of t> feet 
per second. Suppose further that the speed profile within the layer is an 
oscillation of frequency K and amplitude 8. We address the question of how 
this profile can affect the nominal calculation of the ray's horizontal distance 
and transit time through the layer. 

The question is most easily treated using isospeed ray tracing and 
modeling the oscillations as follows: Partition the layer into 2K equithick 
sublayers and assume the sound speed alternates between the values t>+8 and 
d- 8 as we move through these layers. Let 0j be the elevation angle for the 
o+8 speed layers and 02 for the u-8 layers. 

According to isospeed ray tracing formulation, the horizontal distance 
advanced in the layer is 

2 k 

H = ^ E cot K) = Yl cot ^ + cot ^ 2 )l 

and the transit time 

A 2 K -I 

T = ^-Y 

2 K i Vjsinfy) 

_ aJ i i 

2 (u+5)sin(0 1 ) (u-8)sin(0 2 ) 

The nominal values are found by using x> as the speed and 0i as the angle 
throughout the layer. Thus the error in these two values is 
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1 



2 



AH = “[cot(# 2 )-cot( <?,)], 

1 | 

(u + 5)sin(0|) (t>-5)sin(0 2 ) i>sin(0 1 ) 

and 0] and 02 are related by the ray invariant equation 

cos(0i) _ cos(0 2 ) 
u + <5 t>-<5 

The error AH is affected by velocity only through this equation. Notice that 
the errors do not depend upon the frequency K. Using 8/u as the proportion 
of the speed appearing in the amplitude, we can rewrite 

1 1 2 

(l + S/ ujsin^) (1-8/ u)sin(0 2 ) sin^) 





1-8/ v , 1 + 8/ v 2 


2v 

-A 


sin(0j) sin(0 2 ) 

f 1 1 > 


sin(0]) 
(1 + 8/ v) 


2v 


,sin(0 2 ) sin(0!)J 



since the proportion 8/u is believed small. 

Some speculative calculations appear in Table (D-l) for A z = 5 feet and o = 
4800 feet /second. 

The calculation is not very sensitive to values of \), but quite responsive 
to the elevation angle. It should be noted that the signs of AH, AT change if 
the modeled oscillations are in reverse order. This is equivalent to replacing 
with -8. 

This author is not qualified to judge the reality of the suggested values of 
8/ o. Certainly the question deserves more attention. It is common to process 
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some 200 of these 5 foot layers in a ray computation. The error buildup, even 

if the signs change randomly, could be significant, perhaps 15 times AH. 

TABLE D-l. EFFECT OF OSCILLATION 
(v = 4800 feet/second; A z = 5 feet) 



6/u 


0i 


AH 


AT 


.0001 


0.1 


0.48 


-.0001006 




0.15 


0.15 


-.0000301 




0.2 


0.06 


-.0000127 




0.3 


0.02 


-.0000037 




0.4 


0.01 


-.0000015 


.0005 


0.1 


2.18 


-.0004517 




0.15 


0.70 


-.0001432 




0.2 


0.30 


-.0000615 




0.3 


0.009 


-.0000181 




0.4 


0.04 


-.0000074 


.001 


0.1 


3.87 


-.0008032 




0.15 


1.31 


-.0002699 




0.2 


0.58 


-.0001189 




0.3 


0.18 


-.0000357 




0.4 


0.08 


-.0000147 



The support for this calculation proceeds as follows: Suppose X is a 

random variable uniformly distributed on the interval (-AH,AH). Then the 
mean of X is zero and the variance is (AH) 2 /3. If there are n layers to be 
processed then the error in determining the horizontal distance is the sum of 
N independent and identical such X's. It will have mean zero and standard 
deviation AH^n/3; zero plus or minus two standard deviations could be a 
significant amount, especially for small elevation angles. 
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APPENDIX E. CONVERSION OF FOUR TRANSIT TIMES TO INPUTS FOR 



THE RAY TRACING ALGORITHM 



The method in current use for converting the four transit times (ti, t 4 ) 
into an azimuth angle, <j> c , an elevation angle, 0 C and t ac an estimated transit 
time from the source to the acoustic center is outlined and critiqued below. 

The two angles <t> c and 0 C are generated from a description that assumes a 
constant value u for the speed of sound for all points in the array. Then the 
concept of an "apparent position," (X,Y,Z) relative to the acoustic center for 
the sources in a constant speed median is utilized for purposes of estimating 
the two angles. The apparent position must satisfy the system of equations 



This is a system of four equations in three unknowns. Unless the values of tj, 
..., t 4 are singularly coherent, there are an infinity of solutions for X, Y, Z. 

The current policy is to obtain a unique solution by subtracting the fourth 
equation successively from each of the other three, leaving a three by three 
system remaining. I.e., 



(X + D / 2) 2 + ( Y - D / if + (Z - D / 2) 2 = u 2 f, 2 
(X - D / 2) 2 + (Y + D / 2) 2 + (Z - D / 2) 2 = v 2 t\ 
(X - D / 2) 2 + (Y - D / 2) 2 + (Z + D / 2) 2 = u 2 f| 
(X-D/2) 2 + (Y-D/2) 2 + (Z-D/2) 2 = v 2 tj 



(E.l) 



2XD = v 2 



2 YD = v 2 



2 ZD = v 2 




(E.2) 
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and a unique solution. Other rationales supporting this approach can be 
found in [5]. 

At this point an adjustment is made under the name of a direction cosine 
correction (DCC). The quantity D /2 is added to each component of the 
apparent position (in effect making the position relative to the C-phone) and 
the length of the new vector is computed and divided by \)t 4 . Call this 

DCC = ^(X + D/2) 2 +(Y + D/2 ) 2 + (Z + D/2 ) 2 / vt 4 (£.3) 



The translated values are normalized by DCC prior to returning the origin to 
the acoustic center: 



X c = 

v; = 

v c = 



X + D/2 

DCC 

Y + D/2 

DCC 

Y + D/2 

DCC 



-D/2 

-D/2 

-D/2 



(£.4) 



Next come the tilt corrections, required because the corrected apparent 
position is in the coordinate system cs(ac), see Appendix B. This is 
accomplished by applying the approximation of eq. (B.4), i.e., that coefficient 
matrix (instead of the exact one (B.5)) to the vector OC, Y c , Z c ): The Z rotation 
can also be applied at this point, i.e., multiplying by P 3 . Let us call the result of 
all this (X,Y,Z) and form the functions 

sin( 0 c ) = Z / Vx 2 +Y 2 +Z 2 

sin (0 C ) = Y/Vx 2 + Y 2 (£.5) 

cos(</> c ) = X / Vx 2 + Y 2 



and, from these, the angles 0 C and <J) C can be found. 
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Finally, there is need to construct a transit time to the acoustic center 
because it is not measured. Calling the values t ac , the proportionality 
adjustment with t4 is used. 

= f 4 Vx 2 + Y 2 + Z 2 / i/(X + D/2 ) 2 + (Y + D/2) 2 +(Z + D/2 ) 2 (E. 6 ) 
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APPENDIX F. ALTERNATIVE INITIALIZATION OF RAY TRACING 



Current methodology utilizes the four measured ray transit times (tj, t 2 , 
t 3 , t 4 ) from source to the X, Y Z, C hydrophones (respectively), and converts 
them to <t> c , 6 C / and t ac , the azimuth and elevation angles and an estimated 
transit time from source to the acoustic center. It is seen in Section 5 that 
there can be considerable error in this, especially for ()> c . Our alternative is to 
shift the acoustic center to the C-phone and to ignore t 3 - The resulting angles 
will have much smaller error, and the transit time t 4 is used directly. 

We proceed to develop the ray azimuth (spherical coordinate longitude) 
and elevation (spherical coordinate latitude) angles from times ti, t 2 , t 4 , in the 
range coordinate system. (See Appendix B for a general discussion of 
coordinate systems.) The conversion of those times is influenced by the array 
orientation. 

Let Sj,Cj be the sine and cosine of the i ,h Euler angle in the orientation of 
the 3-D sensor array j; i = 1, 2, 3. The conversion of XTILT, YTILT, and ZROT 
into Euler angles is explained in Appendix B, together with methodology for 
locating the positions of the phones at the ends of the array arms. 

Having moved the acoustic center to the C-phone, we are using the 
coordinate system cs(a). Letting Bj represent the j th column of the matrix B, 
eq. (B.3), the locations of the X,Y,Z phones are, respectively, 

D • Bj D • Z?2 D • Bj 
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and, of course, the C-phone is at the origin. These locations in the range 
coordinate system can be found by adding the location of the C-phone as 
indicated in Table B-2. 

Typically there is very little difference in the depths (3 rd component) of 
the X,Y, and C-phones. Hence there is small variability if the speed of sound 
at these three depths and the use of a constant common value, say v, is more 
tenable than the previous use. 

Again we adopt the concept of an apparent position (X,Y,Z) of the sound 
source and the first two direction cosines can be found by applying the law of 
cosines (see the figure) for j = 1, 2. 




where ti, t 2 , t 4 are the signal transit times to the X, Y and C phones, 
respectively. 

For j = 1, the cosine of the marked angle is X/ot 4 and by the law of cosines 
(v ti) 2 = D 2 +(ot 4 ) 2 - 2D • t 4 cos('Fi) 

It follows that 

D 

X= (u • t 4 -o • ti)Co • t 4 + x> • ti)/(2D) 
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and a like equation for j = 2, leading to 



Y = D/2 + Cl)t4 - Ot2> (ot4 + Ut2) / 2D. 

The third direction cosine is obtained from 

cosW = Vl -cos^O-cos^) = -\J 1 - (X2 + Y2) / \) 2 t 



Next we rotate these values into alignment with the test range coordinate 
system, i.e., apply the matrix B, eq. (B.3) 



(F.l) 



The proposed ray elevation (0) and azimuth (<|)p) angles can be found from 

(F.2) 



X' 




'X' 


v c 


= B 


Y 


X. 




Z 



cos(e f ,) = z c / 1 /x c 2 + y c 2 +z ( 



and 

s in(<p p ) = Y c /Jxfi^ 
cos(<ft p ) = X c /Jxf+Y? 

Comment. This technique was designed to treat the speculated shortcomings 
of the procedure currently in use. The results are quite successful, especially 
in reducing azimuth error. It has the curious property of not using 
information provided by the Z-phone. It seems wasteful not to use this 
information. Yet any system that does use it must be required to perform at 
least as well as this one. 
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Appendix G. 
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SUBROUTINE ERRCUR(Hl,Zl,K,A2,XTILT,YTILT,ZROT,D,L,VEL,M, 

* THR,THC,THER,PHR,PHC,PHER,HC,HER,ZC,ZER,TEE,TIMER) 

»»»»»»»»»»» ft.***************************************** »»»»»»»»»»» 

08/10/90 

COMPUTES THE TRUE ELEVATION ANGLES AND THE ESTIMATED ELEVATION 
ANGLES FOR K SOUND SOURCE DISTRIBUTED EQUALLY ON A CIRCLE OF RADIUS 
HI AT DEPTH Z1 FOR A RECIEVING ARRAY AT DEPTH A2 BUT OTHERWISE AT THE 
CENTER OF THE CIRCLE. THE ACCOUSTIC CENTER IS THE GEOMETRIC CENTER OF 
THE ARRAY CUBE. THE CURRENTLY USED METHODOLOGY IS APPLIED. 

»»»»»» *^ »»»»»»» »«^»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»» 



INPUTS: 

HI: RADIUS OF CIRCLE IN FEET 

Z1 : DEPTH OF SOUND SOURCE IN FEET 

K: NUMBER OF SOURCES ON THE CIRCLE 

A2: DEPTH OF THE C-HYDROPHONE 

XTILT/YTILT ,ZROT : ORIENTATION INFORMATION ABOUT THE 
SENSING ARRAY (RADIANS) 

D LENGTH OF ARRAY EDGES. 

L: DEPTH OF LAYER BOUNDARIES. 

M: NUMBER OF RECORDS IN THE VELOCITY DEPTH PROFILE. 

VEL: AVERAGE SPEED OF SOUND IN THE LAYERS. 



OUTPUTS: 

THR: ACTUAL ELEVATION ANGLE (THETA) 

THC: THETA ESTIMATES AT ACCOUSTIC CENTER, CURRENT METHOD 

THER: THETA ERROR AT THE ACCOUSTIC CENTER, CURRENT METHOD 
PHR: ACTUAL AZIMUTH ANGLE (PHI). 

PHG PHI ESTIMATES AT ACCOUSTIC CENTER, CURRENT METHOD. 
PHER: PHI ERROR AT ACCOUSTIC CENTER. 

HO HORIZONTAL ESTIMATE, CURRENT METHOD. 

HER: HORIZONTAL ERROR. 

ZO VERTICAL ESTIMATE, CURRENT METHOD. 

ZER: VERTICAL ERROR. 

T: TRANSIT TIME TO ACCOUSTIC CENTER. 

TIMC: ESTIMATE OF TIME TO ACCOUSTIC CENTER, CURRENT METHOD. 
TIMER: TRANSIT TIME ERROR. 

»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»*»*»»» 



DIMENSION B(5,3),DZ(300),THER(30),HD(5),HC(30),ZC(30) 
DIMENSION L(300),LM(300),LL(300),PX(30),PY(30),T(5),TEE(30) 
DIMENSION THC(30),THR(30),V0(300),V1(300),VEL(300),VV(300) 
DIMENSION HX(5),HY(5),HER(30),ZER(30),TIMER(30) 
DIMENSION X0(3),PHC(30),PHR(30),PHER(30),TIMC(30) 

RE AL*8 A1 ,A2, A2M, ANG,B,C1 ,C2,C3,D,DP,DW,DZ,THER,GG 
REAL*8H,H0,H1,HD,L,LM,LL,LPZ,PIE,P1,P2,PX,PY,S1,S2,S3 
RE AL*8 T,T0,THEC,TH0,TH1,THC,THR,HC,ZC,HER,ZER 
REAL*8 V, V0,V1 , VE L, VV,X,XO,XTILT,Y,YTILT ,Z,Z0,Z1 ,ZROT 
REAL*8CC1,CC2,VV1,VV0,HX,HY,SC,DCC,RAC,RC,TEE 
RE AL*8 C AZ,S AZ,PHC,PHR,PHER,LH,C,T3P,TIMC,TI MER 
REAL*8 SR,SRER,LB,SV,SVU2,SVU,SU2,SU4,G1,G 

PIE = 3.141 59265359D0 
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IEST = 0 
Ml = M 



C DISTRIBUTE THE K SOURCES EQUALLY AROUND THE CIRCLE COUNTER 
C CLOCKWISE FROM THE EAST. 

DO 10I = 1,K 

ANC = 2*PIE*(I-1)/K 
PX(I) = Hl*DCOS(ANG) 

PY(I) = H1*DSIN(ANG) 

PHR(I) = ANG 

IF(ANG.GT.PIE) PHR(I) = ANG - 2*PIE 
10 CONTINUE 

C FORM SINES AND COSINES OF ALL THE EULER ANGLES:ROLL,PITCH,YAW 

52 = DSIN(XTILT) 

C2 = DSQRT(1 - S2**2) 

SI = DSIN(YTILT)/C2 
Cl = DSQRT(1 - SI ”2) 

53 = -DSIN(ZROT) 

C3 = DCOS(ZROT) 

C IN THE COORDINATE SYSTEM HAVING CENTER AT THE C-HYDROPHONE 
C AND POSITIVE-UPWARD, THE LOCATIONS OF THE FOUR HYDROPHONES 
C (RELATIVE TO THE ARM LENGTH D) ARE DEVELOPED NEXT. THIS IS 
C THE TRANSPOSE OF THE MATRIX B IN APPENDIX B. 

B(l,l) = C2*C3 
B(l,2) = C2*S3 
B(l,3) = S2 

B(2,l) = -S1*S2*C3 - Cl»S3 
B(2,2) = -S1*S2*S3 + Cl*C3 
B(2,3) = Sl*C2 
B(3,l ) = -C1*S2*C3 + S1*S3 
B(3,2) = -Cl*S2*S3 - Sl*C3 
B(3,3) = C1*C2 

C LIKE NOTATION WILL BE USED TO LOCATE THE C-HYROPHONE AND THE 
C ACCOUSTIC CENTER (AC). 

DO 12 J = 13 
B(4,J) = 0.0D0 

B(5,J) = 0.5*(B(1,J> + B(2J) + B(3,])) 

12 CONTINUE 



A1 = 0.0D0 
P2 = Z1 



C LOCATE THE HYDROPHONE HORIZONTAL COMPONENTS IN THE COORDINATE 
C SYSTEM CENTERED AT AC. 

DO 14 J = 1,5 

HX(J) = D*(B(J,1) - B(5,l)) 

HY(J) = D*(B(J,2) - B(5,2)) 

14 CONTINUE 

C DETERMINE THE DEPTHS OF THE FOUR HYDROPHONES AND THE AC. 

HD(1) = A2 + D*(B(5,3) - B(l,3)) 

HD(2) = A2 + D*(B(5,3) - B(2,3)) 

HD(3) = A2 + D*(B(5,3) - B(3,3)) 

HD(4) = A2 + D*(B(5,3) - B(4,3)) 

HD(5) = A2 
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C FIND THE DEEPEST HYDROPHONE 
A2M = 0.D0 
DO 51 J=l,4 

IF(HD(J).GT.A2M) A2M = HD(J) 

51 CONTINUE 

C FORM THE SET OF LAYER MIDPOINTS. 

DO 105 1 = 1,M-1 

LM(I) = 5*(L(I) + L(I+1)) 

105 CONTINUE 

LM(M) = LM(M-l) + L(M) - L(M-l) 

C FORM DEPTH INCREMENTS, AND ALL SOUND VELOCITY SLOPES AND 
C INTERCEPTS. 

DO 110 I=1,M-1 

DZ(I)=LM(I+1)-LM(I) 

V0(I)=(LM(I+1)*VEL(I) - LM(I)*VEL(I+1))/DZ(I) 

V1(I)=(VEL(I+1)-VEL(I))/DZ(I) 

110 CONTINUE 
C 

IF(A2M.LT.LM(M-2)) GOTO 126 

C IF A2M IS DEEPER THAN THE LAST LAYER MIDPOINT, THEN WE EXTRAPOLATE 
C THE SOUND VELOCITY PROFILE BY USING A QUADRATIC FUNCTION OVER THE 
C DEEPEST 100 FEET. 

C FIRST COUNT THE NUMBER OF LAYERS (OF THICKNESS DZ(M-2)) TO 
C BE ADJOINED. ALSO MUST EXTEND THE L ARRAY. 

K0 = 2 + MAX(0,NINT((A2M-LM(M-1))/DZ(M-1))) 

MC = 21 

C FIND THE AVERAGE DEPTH FOR THE LAST 100 FEET. 

LB = 0.0D0 

DO 200 J = M+1-MC,M 
200 LB = LB + LM(J) 

LB = LB/MC 

C FORM SUM OF POWERS AND PRODUCTS. 

SV = 0.0D0 
SVU2 = 0.0D0 
SVU = 0.0D0 
SU2 = 0.0D0 
SU4 = 0.0D0 
G1 = 0.0D0 

DO 210 J = M+1-MC,M 
U = LM(J) - LB 
SV = SV + VEL(J) 

SVU = SVU + U*VEL(J) 

SVU2 = SVU2 + U”2 * VEL(J) 

SU2 = SU2 + U**2 
SU4 = SU4 + U**4 
G1 = G1 + V10) 

210 CONTINUE 

G1 = G1/MC 
G = SVU/SU2 

GG = (MC*SVU2 - SU2*S V)/ (SU4*MC - SU2”2) 
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IF(GG.LT.O.ODO) GG = O.ODO 
IF(Vl(M-l).LT.O.ODO) Vl(M-l) = G1 

C PERFORM THE EXTRAPOLATION. 

DO 125 I=M,M+K0 

V1(I) = Vl(I-l) + GG*DZ(M-1) 

LM(I+1) = LM(I) + DZ(M-l) 

VEL(I+1) = VEL(I) + DZ(M-1)»V1(I) 

V0(I) = (LM(I+1)*VEL(I) - LM(I)*VEL(I+1))/DZ(M-1) 

L(I+1) = L(I) + DZ(M-l) 

DZ(I) = DZ(M-l) 

125 CONTINUE 

C UPDATE M, THE NUMBER OF LAYERS 
M = M+KO 

126 CONTINUE 

C THE OUTER LOOP WILL PERFORM COMPUTATIONS FOR THE K SOUND 
C SOURCES. 

C ADJUST DV TABLE TO 25 FT. INCREMENTS. 

CALL VELMOD(L,VEL,Ml,LL,VV,MM) 

C RAYFITTING 

C THE INNER LOOP WILL FIT RAYS TO THE FOUR HYDROPHONES 
C AND THE AC IN THE ORDER X,Y,Z,C AND AC. 

DO50I = l,K 

WRITE!*,*)’ OUTER LOOP I = ’,1/ K = ’,K 
D035 J = 1,5 

PI = DSQRT((PX(I)-HX(J))**2 + (PY(I)-HY(J))**2) 

ZO = HD(J) 

CALL R A YFIT1 ( A1 ,Z0,P1 ,P2,M, V EL,LM,DZ,V 0, VI ,T0,TH0, 

* TH1,IEST) 

C COLLECT THE FIVE TRANSIT TIMES. 

T(J) = TO 

C IN THIS PROGRAM WE KEEP ONLY THE TRUE ELEVATION ANGLE AT AC. 

THR(I) = THO 
35 CONTINUE 
C INNER LOOP COMPLETED. 

C LOCATE THE WATER LAYER, N, CONTAINING THE ARRAY. USE THIS 
C TO DEVELOP THC, THE CURRENTLY USED ESTIMATE OF THO. 

N = MM 
DO 37 J = 2,MM 

IF((LL(J-1).LE.A2).AND.(LL(J).GT.A2)) N = J-l 
37 CONTINUE 

C 

V = VV(N) 

C USE THE FOUR TRANSIT TIMES TO PRODUCE ESTIMATES OF THE 
C ENTRANCE ANGLE. CALCULATE THE PRE-TILT CORRECTED APPARENT 
C POSITION AND INCLUDE THE DIRECTION COSINE CORRECTION. 

DO40J = U 

X0(J)=((V*T(4)-V*T(J))*(V*T(4)+V*T(J)))/(2*D) 

X0(J) = 0.5*D + X0(J) 

40 CONTINUE 
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SC = V*T(4) 

DCC = (DSQRT(X0(1)**2 + X0(2)**2 + X0(3)**2))/SC 
D041 J = 13 

X(XJ) = (XOO)/DCC) - 0.5'D 
41 CONTINUE 

C NEXT MAKE TILT CORRECTIONS 
X = X0(1) - X(K3)*DSIN(XTILT) 

Y = X0(2) - XCX3)*DSIN(YTILT) 

Z = X0(3) + X0(1 )*DSIN(XTILT) + X0(2)»DSIN(YTILT) 

RAC = DSQRT(X**2 + Y**2 + Z**2) 

RC = DSQRT((X + D/2)**2 + (Y + D/2)»*2 + (Z + D/2)**2) 

TO = T(4)*RAC/RC 
TIMC(I) = TO 

C PERFORM Z ROTATION IN THE (X,Y) PLANE. 

X0(1) = C3*X - S3*Y 
X0(2) = S3*X + C3*Y 

C COMPUTE THEC: THE ESTIMATE OF THETA CURRENTLY IN USE. 
THEC = DASIN(Z/DSQRT(X**2 + Y**2 + Z**2)) 

C COMPUTE SINES AND COSINES OF AZIMUTH. 

SAZ = X0(2)/DSQRT(X**2 + Y**2) 

CAZ = X0(1)/DSQRT(X**2 + Y”2) 

PHC(I) = DATAN2(SAZ,CAZ) 

PHER(I) = PHC(I) - PHR(I) 

IF(ABS(PHER(I)).GT.PIE) PHER(I) = PHC(I) + PHR(I) 

IFG = 0 



C RAYTRACE BY THE ISOSPEED METHOD. 

CALL ISOSPEED( Al^O^mTHEC^L^VV^M^^THl ) 

HC(I) = H 
ZC(I) = Z 

C NOW HNISH THE OUTPUT. 

THC(I) = THEC 

C AND THE ERRORS 
TEE(I) = T(5) 

THER(I) = THC(I) - THR(I) 

TIMER(I) = TIMC(I) - T(5) 

HER(I) = HC(I) - HI 

ZER(I) = ZC(I) - Z1 

SR = DSQRT(H1”2 + (A2 - Zl)**2) 

SRER = DSQRT(HC(I )**2 + (A2 - Zl)**2) - SR 
50 CONTINUE 
C OUTER LOOP COMPLETED! 

RETURN 

100 FORMAT(3(5X,E13.6)) 

120 FORM AT(3(5X,F1 5.12)) 

130 FORMAT(3(F12.8,2X)) 

END 
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****** I*************************************************** ****** 

LIBRARY FILE: LIB12.FOR 2/22/90 

*************************************************************** 

*************************************************************** 

SUBROUTINE ISOGRAD1(A1,A2,TO,THO,N,LM,VEL,VO,V1,DZ,H,Z,TH1) 
£**************************************************************** 

C 09/25/89 

C TO : TRANSIT TIME (SEC). 

C THO : ELEVATION ANGLE AT SENSOR (RAD). 

C A1 : HORIZONTAL COORDINATE OF SENSOR. 

C A2 : VERTICAL COORDINATE OF SENSOR, POSITIVE DOWN. 

C V0,V1 : ARRAYS CONTAINING SOUND VELOCITY PARAMETERS. 

C LM : ARRAY CONTAINING LAYER MIDPOINTS. 

C N : INDEX OF DEEPEST LAYER USED. 

C 

£**************************************************************** 

DIMENSION LM(300),V0(300),V 1 (300),DZ(300),VEL(300) 

REAL*8 T0,H,H 0,Z, A 1 , A2,TH0,TH 1 , VEL 
REAL*8 LM,V 0,V 1 ,DZ,Q1 ,Q2 
REAL*8 V A2,R, VP2,TH,RV,D W,DT,X,T 
INTEGER N,IS 

I = N 
T = 0.0D0 
TH=TH0 
HO = A1 

VA2=V0(I)+V1(I)*A2 
RV=DCOS(TH)/VA2 
Z = A2 

DZ(N) = Z - LM(N) 

50 IF(V1(I).EQ.0.0) THEN 

DW = DZ(I)/DSIN(TH) 

DT = DW/V0(I) 

H = HO + DW*DCOS(TH) 

TH1 = TH 
ELSE 

Q2=-V0(I)/V1(I) 

IF (Q2) 51,52,53 

51 IS = -1 
GOTO 54 

52 IS = 0 
GOTO 54 

53 IS = 1 

54 CONTINUE 

Q1=H0 + (Q2-Z)*DTAN(TH) 

R=DSQRT((Q2-Z)**2 + (Ql-H0)**2) 

THl=DACOS(RV*VEL(I)) 

DT=DLOG((DCOS(TH)/(l+DSIN(TH)))*((l+DSIN(THl))/ 

* DCOS(THl)))/Vl(I) 

H=Q1 - IS*R*DSIN(TH1) 

ENDIF 

T=T+DT 

IF (T.GE.TO) GOTO 60 
Z=LM(I) 

HO = H 
TH=TH1 
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1 = 1-1 

GOTO 50 



60 DT=T0+DT-T 

IF(V1(I).EQ.0.0) THEN 
DW = V0(I)*DT 
DZ(I) = DW*DSIN(TH1) 

H = HO + DWDCOS(THl) 

Z = Z-DZ0) 

ELSE 

X=(EXP(DT*V1 (I)))*(l +DSIN(TH) )/ DCOS(TH) 
THl=DACOS((2*X)/ (1 +X**2)) 

H = Q1 - IS*R*DS1N(TH1) 

Z = Q2 - IS*R*DCOS(THl ) 

ENDIF 

C RESTORE THE END LAYERS. 

DZ(I) = LM(I+1) - LM(I) 

DZ(N) = LM(N+1) - LM(N) 

RETURN 
END 

I**************************************************************** 



SUBROUTINE ISOSPEED( A1 ,A2,T0,TH0,L,VEL,M,H,Z,TH1 ) 
£»»»»»*»»»»»»»»»»»*»»»»»»»»»»»*»»»»»»»»»»»»»»»»»»*»»»»**»»»»»»»»» 
C 08/09/89 

C This is a 2-D ray tracing algorithm that mimics the one in 
C Proceedure 5181. It utilizes the assumption that the speed 
C of sound in water is constant for the entire layer encompassed 
C by the layer boundaries. A fixed ray invariant is used 
C throughout the entire migration. 
£»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»**»*»»»*»»»*»»»»*»»»»»»»»»»» 
C INPUTS: 

C A1,A2 - POSITION OF SENSOR (A2>0 IX) WN) 

C TO - TRANSIT TIME 

C THO - ELEVATION ANGLE (OF THE RAY AT THE SENSOR, 
C ALSO CALLED THE ENTRANCE ANGLE) 

C L - ARRAY CONTAINING LAYER BOUNDARIES 
C VEL - ARRAY CONTAINING SOUND VELOCITY AT THE 
C THE LAYER MIDPOINTS 

C OUTPUTS: 

C H,Z - POSITION OF TARGET (SOUND SOURCE) 

C TH1 - ELEVATION ANGLE AT TARGET 
£»»»»*»*»»»»»»»»*»***»»»»»»*»*»»» »»»»*»»»»* ********************** 
DIMENSION L(300),VEL(300) 

REAL*8 A1,A2,C,S,T,DT,DW,TH / DZ,RV / TH1,Z,H,L,VEL,TH0,T0 
Z = A2 
H = A1 
T = 0.0 



C CHOOSE N SUCH THAT L(N) <= A2 < L(N+1). IF A2 IS DEEPER 
C THAN LOWEST LAYER THEN N = M, THE INDEX OF THE DEEPEST LAYER 
C BOUNDARY 
N = M 
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D05I = 2,M 

IF ((L(I-1).LE.A2).AND.(L(I).GT.A2)) N = I - 1 
5 CONTINUE 

RV = DCOS(THO)/VEL(N) 

J = N 
TH = THO 
S = DSIN(THO) 

C = DSQRT(1 - S**2) 

10 DZ = Z - L(J) 

C COMPUTE THE INCREMENTAL SLANT RANGE 
DW = DZ/S 

C COMPUTE THE INCREMENTAL TRAVEL TIME 
DT = DW/VEL(J) 

C ACCUMULATE TOTAL TRAVEL TIME AND TEST 
T = T + DT 

IF (T.GE.TO) GOTO 50 

C UPDATE THE HORIZONTAL AND VERTICAL ACCUMULATIONS 
H = H + DW*C 
Z = Z - DZ 



C USE SNELL'S LAW TO UPDATE THE LAYER ENTRANCE ANGLE AND 
C THE TRIG FUNCTIONS. 

J = J-1 

C = RV * VEL(J) 

S = DSQRT(1 - C**2) 

GOTO 10 



50 T = T - DT 
DT = TO - T 
DW = VEL(J)*DT 
DZ = DW*S 
H = H + DW*C 
Z = Z - DZ 
TH1 = DASIN(S) 

RETURN 

END 

**************************************************************** 



SUBROUTINE RAYFIT1(A1,A2,P1 / P2 / M / VHL,LM,DZ,V0 / VLT0 / TH0, 
TH1,IEST) 

************************************************************** 

09/12/89 

NEW SUBROUTINE TO REPLACE TGEN, RAYTRACING ALGORITHM. 
************************************************************** 



INPUTS: 

A1,A2 - POSITION OF SENSOR (A2 > 0 DOWN) 

P1,P2 - POSITION OF SOUND SOURCE ( P2 > 0 DOWN ) 

LM - ARRAY CONTAINING LAYER MIDPOINTS 
M - NUMBER OF LAYER MIDPOINTS 
VEL - ARRAY CONTAINING SOUND VELOCITY AT THE 
LAYER MIDPOINTS. 
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C VO - SPEED INTERCEPT VALUES 

C VI - SPEED SLOPE VALUES 

C DZ - DEPTH INCREMENTS 

C IEST - FLAG FOR INITIALIZING THE ANGLE 
C OUTPUTS: 

C TO - TRANSIT TIME 

C THO - ELEVATION ANGLE AT THE SENSOR 

C TH1 - ELEVATION ANGLE AT THE SOUND SOURCE 

Z"' »»»»»»»»»»»»»»»»»»»»»»» 



DOUBLE PRECISION VEL(300),DZ(300),LM(300),V0(300) 
DOUBLE PRECISION V1(300),ANG(300),G(300) 

RE AL*8 A1 , A2,P1,P2,T0,TH0,TH1 ,EP,S,C 
REAL*8 H,H0,DW,V A2,VP2,GG,R,Z,TH,RV,Q1 ,Q2 
INTEGER M,IS 



EP = ID-6 



C DETERMINE LAYERS INVOLVED IN RAY FITTING 
N = M 
J = M 

DO 30 I=1,M - 1 

IF ((LM(I).LE.A2).AND.(LM(I+1).GT.A2)) N=I 
IF ((LM(I).LE.P2).AND.(LM(I+1).GT.P2)) J=I 
30 CONTINUE 

C MAKE END CORRECTIONS FOR THE LAYERS 
DZ(N) = A2 - LM(N) 

DZ(J) = LM(J+1) - P2 

C COMPUTE SPEED OF SOUND AT A2 AND P2 
VA2 = V0(N) + V1(N)*A2 
VP2 = V0(J) + V1(J)»P2 
IF(IEST.NE.O) GOTO 50 

C INITIALIZE THE ELEVATION ANGLE AT THE SENSOR, THO, BY 
C FITTING A STRAIGHT LINE SPEED PROFILE BETWEEN P2 AND A2. 
IF(VEL(N).EQ.VEL(J)) THEN 

THO = DATAN((A2-P2)/(P1-A1 )) 

ELSE 

Q2 = (VEL(N)’LM(J) - VEL(J)*LM(N)) / (VEL(N)-VEL(J)) 

Ql = 0.5*(P1+A1)+(0.5*(P2-A2)*(P2+A2-2*Q2))/(P1-A1) 

THO = DAT AN((Q1-A1 )/(Q2- A2)) 

ENDIF 

C OUTER LOOP: SET UP RAY FITTING FOR THO = ELEVATION ANGLE 
50 S = DSIN(THO) 

C = DSQRTU.O - S**2) 

I = N 

RV = C/VA2 
HO = A1 
Z = A2 

60 IF(V1(I).EQ.0.0) THEN 
DW = DZ(I)/S 
H = HO + DW*C 
ELSE 

Q2 = -V0(I)/V1(I) 
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IF (Q2) 61,62,63 
IS = -1 
GOTO 64 

62 IS = 0 
GOTO 64 

63 IS = 1 

64 CONTINUE 

Q1 = HO + (Q2-Z)*S/C 
R = DSQRT((Q2-Z)**2 + (Q1-H0P*2) 

C = RV*VEL(I) 

S = DSQRT(1.0-C**2) 

H = Ql - IS*R*S 
ENDIF 

IF (I.EQ.J) GOTO 80 
H0 = H 
Z = LM(I) 

1 = 1-1 
GOTO 60 

80 TH1 = DACOS(RV*VP2) 

C FRACTIONAL LAYER CORRECTION 

IF(V1(J).NE.0.0) H = Ql - IS*R*DSIN(TH1) 

IF (ABS(H-Pl).LT.EP) GOTO 90 

C RE-ESTIMATE THO 

THO = DATAN(DTAN(TH0)*H/P1) 

GOTO 50 

C PREPARE FOR COMPUTATION OF TRANSIT TIME. 

C COLLECT EXIT AND ENTRANCE ANGLES. 

90 ANG(J) = TH1 
ANG(N+1) = THO 
DO 95 I = J+1,N 

ANG(I) = DACOS(RV*VEL(I)) 

95 CONTINUE 

C COMPUTE TRANSIT TIME 
TO = 0.0D0 
DO 100 1 = J,N 

IF(Vld).EQ.O.O) THEN 

TO = TO + DZ(I)/(V0(I)*DSIN(ANG(I))) 

ELSE 

T0=T0 + DLOG((DCOS( ANG0+1 ))*(1 +DSIN( ANG(I))))/ 

* ((l+DSIN(ANG(I+l))rDCOS(ANG(I))))/Vl(I) 

ENDIF 

100 CONTINUE 

C REMOVE THE END CORRECTIONS. 

DZ(J) = LM(J+1) - LM(J) 

DZ(N) = LM(N+1) - LM(N) 

IEST = 1 

RETURN 
END 

****** *** ******************************************************* 
**************************************************************** 
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SUBROUTINE VELMOD(L,VEL,M,LL,VV,MM) 
£»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»*»»»»»»»»»»»»»»»»»»»»»»»»» 

C 02/22/90 

C THIS PROGRAM TAKES THE VELOCITY SOUND PROFILE GIVEN IN FIVE (5) 

C FOOT INCREMENTS AND CONVERTS IT INTO TWENTY FIVE (25) FOOT INCREMENT 
C PROFILE. 

£»»»»»*»»*»*»*»»»»»»»****»»»»»»*»»*»»»»»»»»»»»»»**»»»»»»»»*»»»»»»*»»»» 

C 

C INPUT: 

C L: DEPTH IN 5 FT INCREMENTS 

C VEL: SOUND VELOCITY IN 5 FT. INCREMENTS 

C M: NUMBER OF ELEMENTS IN DEPTH ARRAY 

C OUTPUT 

C LL: DEPTH IN 25 FT INCREMENTS 

C VV: SOUND VELOCITY IN 25 FT INCREMENTS 

C MM: NUMBER OF ELEMENTS IN DEPTH ARRAY 

£»»»»»»»»»»»»»»»»*»»»»»»»»»»»»»»»»»»»» »*»» it*************************** 



DIMENSION L(300),LL(300),VEL(300),VV(300) 

REAL*8 L,LL,VEL,VV,VS 

MM = INT(M/5) 

MREM = M - 5*MM 
VS = 0.0D0 
DO 10 J = 1,MM 
LL(J) = L(5*J-4) 

W(J) = 0.2*(VEL(5*J-4) + VEL(5*J-3) + VEL(5*J-2) + 
* VEL(5*J-1) + VEL(5*J)) 

10 CONTINUE 

DO 20 J = 1,MREM 

VS = VS + VEL(M+1-J) 

20 CONTINUE 

LL(MM+1) = L(5*MM + 1) 

VV(MM+1) = VS/MREM 

RETURN 

END 



SUBROUTINE ERRPROP(Hl,Zl,K,A2,XTILT / YTILT,ZROT,D,L / VEL,M / 

* THR,THCER,TH1ER,PHR,PHER / PH1ER / H2ER / Z2ER,T4) 
£*»»»*»»****»*»*»*»»**»»»»»»»*******»»*»*»»»»»»»»»*»»»»»»»»*»»»»» 

C 08/22/90 

C POSITION ERROR ANALYSIS WHEN THE ORIGIN IS OVER THE C-HYDROPHONE 
C AND THE PROPOSED SYSTEM IS APPLIED. COMPUTES THE TRUE ELEVATION 
C ANGLES AND ESTIMATED ELEVATION ANGLES FOR K SOUND SOURCES 
C DISTRIBUTED EQUALLY ON A CIRCLE OF RADIUS HI AT DEPTH Z1 FOR A 
C RECIEVING ARRAY AT DEPTH A2 BUT OTHERWISE AT THE CENTER OF THE 
C CIRCLE. A2 IS DEPTH OF THE GEOMETRIC CENTER OF THE ARRAY CUBE. THE 
C ACCOUSTIC CENTER IS THE C-HYDROPHONE. BECAUSE OF DIRECT 
C MEASUREMENT THERE IS NO TRANSIT TIME ERROR. 

C 

£»»»»»»»»»»»»»*»»»»»»»»»»***»»»»»»»»»»»»»»*»»»»»»»»»»»»»»*»*»»»»» 

C 

C INPUTS: 

C HI: RADIUS OF CIRCLE IN FEET 
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C Z1 : DEPTH OF SOUND SOURCE IN FEET 

C K: NUMBER OF SOURCES ON THE CIRCLE 

C A2: DEPTH OF THE CENTER OF THE ARRAY 

C XTILT,YTILT,ZROT: ORIENTATION INFORMATION ABOUT THE 
C SENSING ARRAY (RADIANS) 

C D LENGTH OF ARRAY EDGES. 

C L DEPTH OF LAYER BOUNDARIES. 

C M: NUMBER OF RECORDS IN THE VELOCITY DEPTH PROFILE. 

C VEL: AVERAGE SPEED OF SOUND IN THE LAYERS. 

C 

C OUTPUTS: 

C THR: ACTUAL ELEVATION ANGLE (THETA) 

C THONE: THETA ESTIMATE, PROPOSED METHOD 

C THER: THETA ERROR, PROPOSED METHOD 
C PHR: ACTUAL AZIMUTH ANGLE (PHI). 

C PHI: PHI ESTIMATES, PROPOSED METHOD. 

C PH1ER: PHI ERROR, PROPOSED METHOD. 

C HC: HORIZONTAL ESTIMATE, CURRENT METHOD. 

C HER: HORIZONTAL ERROR. 

C ZC VERTICAL ESTIMATE, CURRENT METHOD. 

C ZER: VERTICAL ERROR. 

C T: TRANSIT TIME TO THE C-HYDROPHONE. 



-'**************************************************************** 



DIMENSION B(5,3),DZ(300),HD(5),HX(5),HY(5),ITH2(20),L(300) 
DIMENSION LM(300),PH1 (20),PH1ER(20),PHC(20),PHER(20) 
DIMENSION PHR(20),PX(20),PY(20),T(4),TH1 ER(20),TH2(20) 
DIMENSION TH2ER(20),THC(20),THCER(20),THONE(20),THR(20) 
DIMENSION V0(300), VI (300) ,VEL(300),X0(30),PZ(2,2),PC(2) 
DIMENSION PP(2,3),E(3),TH3(20),TH3ER(20),TT3(3) 

DIMENSION IFL(2),SSC(2),SP(2),H2(20),Z2(20),H2ER(20) 

DIMENSION Z2ER(20),T4(30) 

REAL*8 B,DZ,HD,HX,HY,L,LM,PH1,PH1ER,PHC,PHER,PHR,PX 
RE AL*8 PY,T,TH1ER,TH2,TH2ER,THC,THCER,THONE,THR,VO,V1 
REAL*8 VEL,XO,PONE,T4 

REAL*8 A1,A2,A2M,ALPH1,ANG,C,C1,C2,C3,CAZ,CT,CX,CX0,CY,CY0 
REAL*8 CZ,CZ0,D,DR,DR1,DP1,DP2,DTDS,EPS,F,FT,GG,H1,HH1,P1 
REAL*8 P2 ,PIE,Q1 ,Q2,Q1 P,R0,R1 ,S,S0,S 1 ,S2,S3,S AZ,ST,SS,T0,T3 
REAL*8 T3P,TH0,TH 1 ,THEC, V, VV0,VV1 ,X,XTILT,Y,Y 0,Y1 ,YTILT 
REAL*8 Z,Z0,Z1 ,ZROT,ZZl ,SC,DT 

REAL*8 PP,E,DEL,A,TH3,TH3ER,TT3,SSC,SP,PZ,PC 
REAL*8 A11,A12,A21,A22,B1,B2,H2,H2ER,Z2,Z2ER,THE 
REAL*8 LB,SV,SVU,SU2,SU4,SVU2,G,U 
INTEGER ITH2,IFL 

PIE = 3.14159265359D0 
EPS = ID-6 
IEST = 0 
Ml = M 



C DISTRIBUTE THE K SOURCES EQUALLY AROUND THE CIRCLE COUNTER- 
C CLOCKWISE FROM THE EAST. 

DO 101 = 1,K 

ANG = 2*PIE*(I-1)/K 



88 



PX(I) = HrDCOS(ANC) 

PY(I) = H1*DSIN(ANG) 

PHR(I) = ANG 

IF(ANG.GT.PIE) PHR(I) = ANG - 2*PIE 
10 CONTINUE 

C FORM SINES AND COSINES OF ALL THE EULER ANGLES:ROLL,PITCH,YAW 

52 = DSIN(XTILT) 

C2 = DSQRT(1 - S2**2) 

SI = DSIN(YTILT)/ C2 
Cl = DSQRT(1 - Sl**2) 

53 = -DSIN(ZROT) 

C3 = DCOS(ZROT) 

C IN THE COORDINATE SYSTEM HAVING CENTER AT THE C-HYDROPHONE 
C AND POSITIVE-UPWARD, THE LOCATIONS OF THE FOUR HYDROPHONES 
C (RELATIVE TO THE ARM LENGTH D) ARE DEVELOPED NEXT. 

B(l,l) = C2‘C3 
B(l,2) = C2*S3 
B(13) = S2 

B(2,l) = -SrS2*C3 - C1*S3 
B(2,2) = -S1*S2*S3 + C1*C3 
B(2,3) = S1*C2 
B(3,l) = -C1*S2*C3 + S1*S3 
B(3,2) = -C1*S2*S3 - S1*C3 
B(3,3) = CUC2 

C LIKE NOTATION WILL BE USED TO LOCATE THE C-HYROPHONE AND THE 
C ARRAY CENTER. 

DO 12 J = 13 
B(4,J) = 0.0D0 

B(5J) = 0.5*(B(J,1) + B(],2) + B(J3)) 

12 CONTINUE 



A1 = 0.0D0 
P2 = Z1 



C LOCATE THE HYDROPHONE HORIZONTAL COMPONENTS IN THE COORDINATE 
C SYSTEM CENTERED AT C-HYDROPHONE. 

DO 14 J = 13 

HX(J) = D*B(Jd) 

HY(J) = D*B(J,2) 

14 CONTINUE 

C DETERMINE THE DEPTHS OF THE FOUR HYDROPHONES AND THE ARRAY CENTER. 
HDd) = A2 + D*(B(53) - B(l,3)) 

HD(2) = A2 + D»(B(53) - B(2,3)) 

HD(3) = A2 + D*(B(5,3) - B(3,3)) 

HD(4) = A2 + D*(B(5,3) - B(4,3)) 

HD(5) = A2 

C FIND THE DEEPEST HYDROPHONE 
A2M = 0.D0 
DO 51 J=l,4 

IF(HD(J).GT.A2M) A2M = HD(J) 

51 CONTINUE 

C FORM THE SET OF LAYER MIDPOINTS. 
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DO 105 1 = 1,M-1 

LM(I) = .5*(L(1) + L(1+1)) 

105 CONTINUE 

C FORM DEPTH INCREMENTS, AND ALL SOUND VELOCITY SLOPES AND 
C INTERCEPTS. 

DO 110 I=l,M-2 

DZ(I)=LM(I+1)-LM(I) 

V0(I)=(LM(I+1)*VEL(I) - LM(I)*VEL(I+1))/DZ(I) 

V1(I)=(VEL(I+1)-VEL(I))/DZ(I) 

110 CONTINUE 

LM(M) = LM(M-l) + DZ(M-2) 

C 

I F( A2M .LT.LM ( M-l )) GOTO 126 

C IF A2M IS DEEPER THAN THE LAST LAYER MIDPOINT, THEN WE EXTRAPOLATE 
C THE SOUND VELOCITY PROFILE BY USING A QUADRATIC FUNCTION OVER 
C THE DEEPEST 100 FEET. 

C FIRST COUNT THE NUMBER OF LAYERS (OF THICKNESS DZ(M-2)) TO 
C BE ADJOINED. ALSO MUST EXTEND THE L ARRAY. 

K0 = 2 + MAX(0,NINT((A2-LM(M-l))/DZ(M-2))) 

C FIND AVERAGE DEPTH OF LAST 100 FEET. 

LB = 0.0D0 
DO 43 I = M-21,M-1 
LB = LB + LM(I) 

43 CONTINUE 
LB = LB/21 

C FORM SUMS OF POWERS AND PRODUCTS. 

SV = 0.0D0 
SVU2 = 0.0D0 
SVU = 0.0D0 
SU2 = 0.0D0 
SU4 = 0.0D0 
DO 45 I = M-21,M-1 
U = LM(I) - LB 
SV = SV + VEL(I) 

SVU = SVU + U*VEL(I) 

SVU2 = SVU2 + U**2 * VEL(I) 

SU2 = SU2 + U**2 
SU4 = SU4 + U**4 
45 CONTINUE 

G = SVU/SU2 

GG = (2.1-SVU2 - SU2*SV)/(SU4 - SU2”2) 

Vl(M-l) = G 

C PERFORM THE EXTRAPOLATION. 

DO 125 I=M,M+K0 

Vl(I-l) = Vl(I-2) + GG*DZ(M-1) 

LM(I) = LM(I-l) + DZ(M-2) 

VEL(I) = VEL(I-l) + DZ(M-2)*V1(I-1) 

V0(I-1) = (LM(I)*VEL(I-1) - LM(I-l)*VEL(I))/DZ(M-2) 

L(I+1) = L(I) + DZ(M-2) 

DZ(I-l) = DZ(M-2) 
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125 CONTINUE 



C UPDATE M, THE NUMBER OF LAYERS 
M = M+KO 
126 CONTINUE 

C LOCATE THE WATER LAYER, N, CONTAINING THE ARRAY. 

N = M 

DO 37 J = 2,M 

IF((LM(J-1).LE.HD(4)).AND.(LM(J).GT.HD(4))) N = J-l 
37 CONTINUE 
C 

V = V0(N) + V1(N)*HD(4) 

C THE OUTER LOOP WILL PERFORM COMPUTATIONS FOR THE K SOUND 
C SOURCES. 

C RAYFITTING 

C THE INNER LOOP WILL FIT RAYS TO THE FOUR HYDROPHONES 
C IN THE ORDER X,Y,Z, AND C. 

DO 50 1 = 1,K 

WRITE!*,*)' OUTER LOOP I = ',1/ K = ’,K 
IW1 =0 
DO 35 J = 1,4 

PI = DSQRT((PX(I)-HX(J))**2 + (PY(I)-HY(J))*»2) 

Z0 = HD(J) 

CALL RA YFIT1 ( A 1,Z0,P 1 ,P2,M,V EL,LM,DZ, V 0,V 1 ,T0,TH0, 

» TH1,IEST) 

C COLLECT THE FOUR TRANSIT TIMES. 

T(J) = TO 

C IN THIS PROGRAM WE KEEP ONLY THE TRUE ELEVATION ANGLE AT THE 
C C-HYDROPHONE. 

THRO) = THO 
T4(I) = T(4) 

35 CONTINUE 

C CALCULATE THE PRE-TILT CORRECTED APPARENT POSITION 
DO40J = 1,3 

X00)=(D*»2 + (V‘T(4)-V*T(J))*(V»T(4)+V*T(J)))/(2*D) 

40 CONTINUE 

C COMPUTE DIRECTION COSINES . 

CXO = X0(1)/(V*T(4)) 

CYO = X0(2)/(V*T(4)) 

CZO = DSQRT(1 - CX0**2 - CY0**2) 

C PERFORM EXACT TILT CORRECTIONS AND THE ROTATIONAL ALIGNMENT. 
42 CX = B(1,1)*CX0 + B(2,l )*CY0 + B(3,1)*CZ0 

CY = B(1,2)*CX0 + B(2,2)*CY0 + B(3,2)*CZ0 
CZ = B(1,3)*CX0 + B(2,3)*CY0 + B(3,3)*CZ0 
IF0W1.EQ.1) THEN 

TH2(I) = 0.5*PIE - DACOS(CZ) 

GOTO 49 
ENDIF 

SAZ = CY/DSQRT(CX**2 + CY**2) 

CAZ = CX/DSQRT(CX**2 + CY*»2) 
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PH1(I) = DATAN2(SAZ,CAZ) 

PHIER(I) = PH1(I) - PHR(I) 

IF(ABS(PH1ER(I)).GT.PIE) PHIER(I) = PH1(I) + PHR(I) 

THONE(I) = 0.5TIE - DACOS(CZ) 

THIER(I) = THONE(I) - THR(I) 

49 CALL ISOGRAD1(A1 / ZO,TO / THONE(I) / N / LM / VEL / VO,V1,DZ,H2(I) / 

* Z2(I),THE) 

H2ER(I) = H2(I) - HI 
Z2ER(I) = Z2(I) - Z1 

50 CONTINUE 

C OUTER LOOP COMPLETED! 

RETURN 

100 FORMAT(3(5X,E13.6)) 

120 FORMAT(3(5X,F15.12)) 

130 FORMAT(4(F10.8,2X) < F13.8 / 1X,F12.8 / 2X / 2(F12.8,2X),F12.8) 

140 FORMAT(5X,'The transit time to the z-phone is not bracketed') 

END 
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